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ABSTRACT 


A  Computer  program  to  determine  the  ultimate  failure 
load  of  a  circular  bulkhead  with  horizontal  and  vertical 
stiffeners  is  presented.  The  circular  plate  is  treated  as  a 
grillage  with  the  plating  acting  as  an  effective  flange  on 
the  stiffener.  The  grillage  is  composed  of  a  combination  of 
transverse  and  longitudinal  members,  however,  the  spacing 
between  nodal  points  is  not  uniform.  The  plate-st i f f ener 
combination  is  subjected  to  a  uniform  hydrostatic  pressure 
that  is  modelled  as  equal  point-forces  at  each  of  the 
interior  nodes.  A  matrix  system  of  inequality  constraints 
(on  maximum  moment  values)  and  nodal  point  equilibrium 
equations  is  established  and  the  applied  point-force  is 
maximized  using  a  Linear  Programming  maximization  routine 
{ZX3LP)  from  the  International  Mathematical  and  Statistical 
Library  (IMSL).  A  detailed  explanation  of  the  computer 
program  is  presented. 

The  ultimate  failure  load  for  the  unsupported  segments 
of  plating  (between  stiffeners)  is  investigated.  A  brief 
review  of  the  use  of  the  Upper  and  Lower  Bound  Theorems  is 
presented. 

Results  from  the  computer  program  are  compared  to 
experimental  data  for  similarly  stiffened  circular  plating 
subjected  to  a  uniform  pressure  until  collapse. 

Thesis  Supervisor:  Dr.  Joao  Manuel  Gomes  de  Oliveira 
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CHAPTER  1   INTRODUCTION 


BACKGROUND 


In  the  last  decade,  the  demand  for  new  and  more 
efficient  modes  of  marine  transportation  has  forced 
naval  architects  and  marine  engineers  to  search  for 
practical  and  reliable  methods  of  structural  analysis 
which  would  take  full  advantage  of  the  inherent 
strength  of  the  materials  of  construction  in  providing 
a  safe  design  with  a  minimal  structural  weight. 

Elastic  design  methods  provide  an  accurate  and 
reliable  technique  for  predicting  structural 
deformations  resulting  from  applied  design  loadings 
whenever  the  material  behavior  is  linear  elastic. 
However,  in  applications  where  the  total  deformation  is 
not  a  limiting  factor,  elastic  analysis  methods  do  not 
provide  an  effective  technique  for  prediction  of  the 
ultimate  failure  load  of  a  given  structural  design. 
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In  an  attempt  to  take  full  advantage  of  the 
inherent  ductility  of  structural  steels,  a  design 
methodology  was  developed  which  treats  load  levels 
exceeding  that  necessary  to  reach  initial  yield  stress 
(elastic  limit).  This  methodology,  called  plastic 
limit  analysis,  is  not  useful  for  predicting  local 
stress  levels  or  deformations,  but  is  a  reliable  and 
effective  method  of  determining  the  ultimate  collapse 
load.  In  those  cases  where  structural  deformation  is 
not  a  limiting  factor,  or  where  the  ultimate 
load-bearing  capacity  of  a  given  structure  is  a  primary 
concern,  plastic  limit  methods  of  analysis  provide  the 
structural  engineer  with  an  effective  tool  for 
accurately  determining  the  load  limit  of  a  given 
structure.  If  the  designer  requires  additional 
information,  such  as  local  stresses  and  deformations, 
or  if  material  strain  hardening  or  finite  deflections 
are  an  important  consideration,  then  a  numerical 
analysis  using  an  incremental  theory  of  elasticity  will 
need  to  be  performed  {1}. 

A  detailed  review  of  available  literature 
concerning  plastic  limit  analysis  of  structures  with 
stiffened  and  unstiffened  plating   indicated  a  great 
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deal  of  recent  interest  in  applying  "plastic" 
techniques  to  structural  analysis  and  design,  with  a 
view  toward  structural  optimization.  By  far,  the 
greatest  majority  of  st if f ened-plat ing  applications 
developed  utilizing  plastic  limit  analysis  dealt  with  a 
class  of  problems  that  could  most  generally  be 
categorized  as  stiffened  rectangular  plating  with 
uniformly  spaced  stiffeners. 

An  unstiffened  circular  plate  has  been  analyzed 
using  plastic  limit  analysis  {2}  with  a  linear 
programming  technique  and  the  "Lower  Bound"  Theorem  (to 
be  discussed  in  detail).  Similarly,  a  stiffened 
circular  plate  was  studied  using  elastic  analysis  to 
optimize  the  plate-stiff ener  combination  {3}.  In  this 
analysis,  however,  the  stiffeners  were  a  radial  and 
circumferential  system  of  stiffeners  (axisymmetr ic ) 
which  were  compatible  with  the  axisymmetric  circular 
plate,  and  the  stiffened  circular  plate  was  optimized 
by  maximizing  the  load  per  unit  deflection. 

In  the  field  of  submersible  vehicle  design,  a 
cylindrical  hull  form  provides  the  best  compromise 
between   competing   criteria  of  adequate  hull  strength, 
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maximum  internal  volume,  and  minimum  hydrodynamic  drag. 
As  a  result,  submersible  vehicle  bulkheads  are 
essentially  stiffened  circular  plates.  Due  to  the 
requirement  for  minimizing  fabrication  cost,  the 
bulkhead  stiffening  system  is  composed  of  a  grid 
network  of  horizontal  and  vertical  stiffeners,  not 
radial  and  circumferential  stiffeners.  Thus,  the 
as-constructed  submersible  vehicle  bulkhead  cannot  be 
analyzed  by  an  adaptation  of  available  elastic  (or 
plastic)  axisymmetr ic ,  orthotropic  plate  analysis 
methods.  Also,  ultimate  failure  is  the  design 
criteria,  without  any  limit  on  deformation. 

In  their  paper] {4},  Palermo  and  Bart  applied  the 
concept  of  plastic  analysis  to  the  as-constructed 
submersible  vehicle  bulkhead.  The  method  detailed  in 
this  paper  combined  the  concept  of  a  plastic  hinge  and 
elastic  grillage  analysis  to  determine  the  required 
scantlings  for  the  main  horizontal  girder.  However, 
the  vertical  stiffeners  were  analyzed  individually  with 
the  main  horizontal  girder  assumed  to  be  a  fixed 
support  for  each  half  of  the  vertical  stiffener.  More 
recent  methods  allow  for  plastic  analysis  of  the 
grillage  as  a  whole. 
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This  thesis  will  attempt  to  analyze  the  stiffened 
circular  bulkhead  as  a  grillage,  adapting  plastic 
grillage  limit  analysis  techniques  developed  by  Hodge 
{5}  for  a  rectangular  grillage.  The  net  result  will  be 
an  interactive  computer  program  (CIRCPLAT)  designed  to 
be  inexpensive  to  run  and  simple  to  use,  providing  a 
reasonably  accurate,  yet  conservative  estimate  of  the 
collapse  load  of  a  user-defined,  stiffened,  circular 
plate. 

The  results  of  the  CIRCPLAT  analysis  will  be 
compared  to  actual  bulkhead  failure  test  data. 
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CHAPTER  2    PROBLEM  DEFINITION 


FORMULATION  OF  THE  PROBLEM 


The  basic  problem  under  consideration  is  to 
determine  the  failure  load  of  a  built-in,  circular, 
stiffened  bulkhead  using  horizontal  and  vertical 
stiffeners  when  it  is  subjected  to  a  uniform,  lateral, 
hydrostatic  pressure.  To  find  an  accurate  solution  in 
the  simplest,  most  efficient  manner,  it  is  first 
necessary  to  decide  what  is  the  most  effective  way  to 
model  the  problem.  Reference  {6}  provided  a  recent, 
well  researched  compendium  of  current  methods  of  gross 
panel  analyses  along  with  the  inherent  assumptions  and 
limitations  of  each  method,  and  a  comparison  of  their 
relative  accuracy.  Figure  1  provides  a  brief  summary 
of  the  results  of  comparison  of  analyses  methods  given 
by  reference  {6}.  Reference  {6}  also  indicated 
excellent  results  using  plastic  limit  analysis.  Based 
on  a  very  favorable  comparison  between  the  finite 
element   program   results   and  the  "grillage  analysis". 
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and,  the  close  correlation  between  both  methods  and 
experimental  data,  it  appeared  to  be  highly 
advantageous  to  model  the  stiffened  circular  plate  as  a 
grillage  of  intersecting  beams  with  an  effective  width 
of  plating  acting  as  a  flange  to  the  stiffener.  With 
this  model,  it  should  be  possible  to  develop  a  computer 
program  that  is  relatively  inexpensive  to  run,  yet 
sacrifices  little  or  nothing  in  accuracy  for 
determination  of  the  failure  load. 
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Total  Stress  (psi)     Percentage 

Centerline  Difference 

Longitudinal 

Stiffener  at  Fixed        ** 

Boundary  End 

Orthotropic    1st  Order       23,320  -11% 

Plate 

Analysis      2nd  Order       23,585  -10% 

Beam  &  Plate     26,285  +.04% 

Finite        Elements 
Element 
SOLID  SAP     Orthotropic      22,933  -12% 

Plate 

Elements 

Finite  Element  26,186  

ICES  STRUDL 

Grillage  Analysis             25,721  -2% 

Simple  Beam  Theory            30,184  +15% 

**  Percentage  difference  with  respect  to  ICES  STRUDL  Results 

Figure  1.    Stiffened  Panel  Analysis  Method  Comparison. 


GRILLAGE  MODEL:   ASSUMPTIONS 


The  use  of  the  grillage  model  in  conjunction  with 
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plastic   limit   analysis  involves  several  approximating 
assumptions: 


1.  The  in-plane  compressive  loading  is  small  compared 
to  the  lateral  loading  and  its  effect  is 
negligible. 

2.  The  influence  of  the  shear  force  on  the  formation 
of  the  plastic  hinge  is  negligible. 

3.  The  torsional  rigidity  of  the  stiffeners  is  small 
and  its  effect  is  negligible. 

4.  The  effect  of  the  plating  can  be  accounted  for  by 
incorporating  an  "effective  breadth"  of  plating  as 
a  flange  on  the  stiffener. 

5.  Plane  sections  remain  plane  in  the  plastic  range. 

6.  Deflections  are  inf intessimal . 

7.  The  material  is  elast ic-perf ectly  plastic; 

8.  Drucker's  Postulate  is  valid. 


Reference] { 5}  discusses  the  validity  of  assuming 
that  the  effect  of  in-plane  stresses  are  negligible. 
Experimental  evidence  indicates  that  if  the  in-plane 
stress  is  less  than  15  percent  of  the  compressive  yield 
stress,  the  effect  of  the  axial  force  can  be  neglected. 
Current  submersible  pressure  hull  design,  utilizing  a 
ring-stiffened  cylinder,  provides  adequate  support 
against   hull   compressive  effects  sufficient  to  ensure 
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that  compressive  axial  stresses  in  the  bulkhead 
stiffeners  will  not  exceed  15  percent  of  yield  stress. 
References  {4}  and  {6}  contain  correction  factors  that 
can  be  applied  to  the  stiffeners'  plastic  section 
moduli  to  account  for  the  effect  of  in-plane  stresses 
above  15  percent  of  yield.  These  correction  factors 
were  not  incorporated  in  the  CIRCPLAT  program.  The 
program  could  easily  be  modified  to  include  these 
corrections  if  large  axial  forces  were  to  be  included 
in  the  analysis. 

The  assumption  that  shear  force  is  negligible  is  a 
reasonable  one  for  long,  slender  beams.  For  short, 
deep  beams,  a  correction  factor  must  be  applied  to  the 
section  modulus.  To  study  the  relative  magnitude  of 
shear  and  bending  stress  in  a  beam,  consider  a 
rectangular,  prismatic  beam  with  a  single  concentrated 
load  as  shown  in  Figure  2:  (Reference  {6}) 
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Figure  2.    Rectangular  Beam  with  Concentrated  Load. 


The  maximum  bending  moment,  which  occurs  at 
midspan  is:  M=PL/4.  The  maximum  bending  stress  will  be 
at  the  same  point: 

Omax  =  Mc/I  =  (3/2) (PL/bh^ ) . 

The  shear  force  is  constant  at  P/2  between  each 
support  and  the  load  P.  The  formula  for  maximum  shear 
stress  in  a  rectangular  beam  is: 

'i:max=(3/4)(P/bh)  . 
Then  ffmax/Zmax  =  2L/h.   Thus,  we  can  see  that  when  L  is 
many   times   (10  or  more)  greater  than  h,  the  magnitude 
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of   the  maximum  shear  stress  as  compared  to  the  maximum 
bending  stress  is  negligible. 

Reference] {7}  discusses  the  practical  significance 
of  torsional  rigidity  in  grillage  analysis,  and 
concludes  that  the  effect  of  torsion  on  the  plastic 
moment  for  I-beams  is  0.1  percent,  and  is  therefore 
negligible. 

The  question  of  incorporating  an  "effective" 
breadth  of  plating  as  a  flange  on  the  stiffener  is  one 
that  appears  to  have  a  number  of  answers  as  to  how  much 
breadth  to  include  as  "effective."  Reference  {4}  used 
a  value  equal  to  the  stiffener  spacing,  up  to  a  maximum 
of  60t,  where  t=plate  thickness.  Reference  {7} 
utilized  a  formula  developed  by  Schade: 

B  =  l.l/(l+2/r^) 
where  B  is  the  fraction  of  stiffener  spacing  that  is 
effective,  and  r  is  the  ratio  of  the  distance  between 
points  of  zero  bending  moment  and  the  plate's  breadth. 
In  reference  {8},  extensive  experimental  research  by 
Kendrick  was  performed  specifically  to  determine  the 
effective  breadth  of  plating  to  be  included  as  a  flange 
to   the   stiffener  when  calculating  the  plastic  section 
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modulus.  The  results  of  this  research  showed  that 
experimental  plastic  moments  agreed  quite  closely  with 
the  theoretic  values  calculated  by  assuming  all  the 
plating  effective  as  a  flange  to  the  stiffener.  For 
this  reason,  the  plastic  section  modulus  calculation 
incoporated  in  the  CIRCPLAT  program  utilizes  the  full 
stiffener  spacing  as  the  "effective"  breadth,  with  no 
maximum  limit  (such  as  60t). 

The  remaining  assumptions  are  the  standard 
assumptions  utilized  in  deriving  the  equations  for 
plastic  limit  analysis. 


PLATE  MODEL  ASSUMPTIONS 


As  an  additional  consideration,  it  is  necessary  to 
consider  the  possibility  of  plate  failure  between 
stiffeners.  Because  this  is  a  circular  plate 
intersected  by  horizontal  and  vertical  stiffeners,  the 
resulting  sections  of  unsupported  plate  created  come  in 
a  variety  of   shapes   (see   figure  3).   Plastic  limit 
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analysis   can   similarly   be   applied   to   these   plate 
sections  to  determine  their  ultimate  failure  load. 


Figure  3.    Stiffened  Circular  Plate  Geometry. 


In   applying   plastic  limit  analysis  to  thin  plates 
it  is  assumed: 

1.  Elastic-perf ectly  plastic  material 

2.  No  in-plane  loads 

3.  Inf initessimal  deflections 
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4.  No  shear  effects 

5.  Johansen's  yield  surface  applies 

Reference  {9}  details  the  use  of  Upper  and  Lower  Bound 
Theorems  of  Plasticity  for  determining  plate  failure 
load  subject  to  the  above  assumptions. 
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CHAPTER  3  STRUCTURAL  ANALYSIS 


BASIC  PRINCIPLES 


Reference  {1}  presents  an  excellent  development  of 
the  theory  of  plastic  analysis  of  structures.  In 
understanding  plastic  failure  analysis,  it  is  first 
necessary  to  understand  the  concept  of  the  plastic 
hinge.  Consider  pure  bending  of  a  prismatic  beam  of 
rectangular  cross-section  subject  to  an  applied  moment. 
Coordinate  axes  and  sign  convention  are  positive  as 
indicated  in  figure  4. 
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Figure  4 


Rectangular  Beam  subject  to  pur'e  bending 


The  bending  moment  is  given  by: 


m=.bJ 


z(lx(z)dZ 


-H 


The  maximum  elastic  bending  moment  will  occur  when 
the  stress  at  the  outer  fiber  reaches  yield.  Then  the 
moment  in  this  case  becomes: 

Me=(4/3) (BH^Oo) 

If  the  applied  moment  is  increased  still  further, 
the  outer  fibers  will  begin  to  deform  plastically,  and 
the   elastic   to   plastic   interface  will  begin  to  move 


Chapter  3  Structural  Analysis 


23 


toward   the   neutral   axis   until,  eventually,  k=0  (see 
figure  5. ) 


— ^ 
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K- 
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Figure  5.    Elastic-Plastic  Stress  Distribution 


At  the  point  where  k=0,  the  beam  can  not  support 
any  further  increase  in  applied  moment,  and  a  plastic 
hinge  is  said  to  have  developed.  This  fully  plastic 
moment,  denoted  by  Mo,  for  a  rectangular  section,  is: 

Mo  =  2  0oBH 

In  the  partially  plastic  range  between  Me  and  Mo, 
the  deformation  is  theoretically  controlled  by  the 
assumption  that  plane  sections  remain  plane. 
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A  yield  hinge,  then,  has  the  following  properties: 

1.  If  the  applied  moment  is  less  than  Mo,  the  yield 
hinge  will  transmit  the  full  moment  without 
allowing  any  rotation. 

2.  If  the  applied  moment  is  equal  to  Mo  in  magnitude, 
the  yield  hinge  will  transmit  the  full  moment,  but 
will  permit  rotations  of  arbitrary  magnitude  in  the 
direction  of  the  applied  moment. 

3.  The  yield  hinge  will  not  transmit  any  moment 
greater  than  Mo  in  magnitude. 

With  the  concept  of  a  yield  hinge  defined,  failure 
of  a  structure  occurs  when  sufficient  yield  hinges  are 
formed  to  allow  inf initessimal  motion  of  all  or  part  of 
the  structure.  The  existance  of  sufficient  yield 
hinges  is  termed  a  "mechanism". 

Several   important  concepts  involved  in  the  idea  of 

a   failure  mechanism  are  exemplified  in  the  case  of  a 

simple,   indeterminate  structure   under   a  distributed 

load.   Consider   the  case  of  a  beam  fixed  at  one  end. 

Chapter  3  Structural  Analysis  25 


simply   supported  at   the   other   end,   and   under   a 
uniformly  distributed  load  P,  as  in  figure  6. 


/I 


i        i        i 


+  M 


L  - 
(a) 


1 


Figure  6.    Beam  with  Distributed  Load 


First,  it  can  be  seen  that  the  problem  is 
indeterminate,  because  there  are  3  unknown  reactions 
and  only  2  equations  of  equilibrium  (horizontal  force 
equilibrium  is  implicit).  From  elastic  analysis,  we 
know  that  the  shape  of  the  moment  distribution  is  that 
of  figure  6.c.  Thus,  we  know  that  a  plastic  hinge  will 
form  at  the  fixed  end,  and  one  will  form  somewhere  near 
the  center.  As  soon  as  the  first  hinge  is  formed,  the 
problem   becomes   determinate   since   we   have  one  less 
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unknown,    leaving    2   unknowns   and   2   equilibrium 
equations. 

For  example,  if  the  first  hinge  forms  at  the  fixed 
end,  then  summing  moments  about  point  A  gives: 

-Mo+RgL-PL  /2  =  0 
Summing  y-directed  forces  gives: 

R  +R„  =  PL 

A   B 

These  two  equations  can  be  solved  for  both  R^  and 
Rg .  Although  the  problem  is  now  determinate,  we  do  not 
yet  have  a  mechanism.  The  plastic  hinge  at  A  would 
allow  rotation,  but  the  simple  support  at  B  still 
prevents  it,  so  a  second  hinge  must  form  somewhere  near 
the  center  of  the  beam.  These  two  hinges  now  form  a 
collapse  mechanism.  Plastic  hinge  rotation  can  proceed 
unimpeded. 


UPPER  BOUND  THEOREM 


Continuing  with  the  example  of  Figure  6,  assume 
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that  the  second  hinge  forms  exactly  at  the  center  of 
the  beam.  The  internal  work  of  deformation  is  the 
product  of  the  plastic  moment  times  the  angle  through 
which  it  rotates.  For  the  case  of  the  beam  in  Figure 
6,  if  the  hinge  at  the  fixed  end  rotates  by  ©,  the 
center  hinge  rotates  by  26.  Work  is  independent  of 
direction  of  rotation.   Then 

Wi=Mo^+MQ(25)=3MQ0. 

External  work  is  the  product  of  applied  force  times 
the   distance   moved.    For   a   distributed   load,   the 

2 

average  distance  moved  is  1/4(L  d) .  If  the  applied 
load  is  large  enough  to  do  sufficient  external  work  to 
equal  the  required  internal  work  of  the  collapse 
mechanism,    then   collapse   will   occur.    Therefore, 

2 

[p/4](L  d)  =  3Md        will      yield      the   value   of    P   required    for 

2 

our  assumed  mechanism.  This  result,  P=12Mq/L  ,  is  the 
required  distributed  load  to  cause  the  "assumed" 
collapse  mechanism  to  occur.  If  the  assumed  collapse 
mechanism  were  the  exact  collapse  mechanism,  the  value 
for  P  would  be  the  exact  value  to  cause  collapse. 
Looking  at  Figure  6c.,  the  moment  distribution  shows 
that  the  maximum  moment  occurs  to  the  right  of  the 
center   of   the   beam.   As  a  result,  the  actual  plastic 
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hinge  will  form  to  the  right  of  the  center.  As  a 
result,  the  collapse  load,  P,  determined  above,  is  not 
the  exact  collapse  load.  The  Upper  Bound  Theorem 
essentially  states  that  the  failure  load  determined  by 
an  assumed  collapse  mechanism  that  is  not  the  exact 
collapse  mechanism  will  be  an  upper  bound  to  the  exact 
failure  load.  Thus,  the  load,  P  =12Mq/L  ,  is  greater 
than  the  actual  collapse  load,  and  is  an  Upper  Bound. 


LOWER  BOUND  THEOREM 


Remembering  that  once  the  first  plastic  hinge 
forms,   the   problem  becomes  determinate,  the  following 

equations  can  be  written: 

2 
M(0)  =       R^L  -  PL      =  -Mo 

2 
M(Z)  =  Rg(L-Z)-  P(L-Z)/2  =   Mo 

M'  (Z)  =    -Rb  +  P(L-Z)    =   0 

where   Z  is  the  point  at  which  the  second  plastic  hinge 

forms.    With   3  equations  and  3  unknowns,  the  reaction 

R,  distance  Z  and  load  P  can  be  determined.   The  result 

is: 
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R   =    4.8284  Mo/L 
P   =  11.65685  Mq/L^ 
Z   =     0.58579  L 

From  the  moment  curve  of  Figure  6c,  it  can  be  seen 
that  if  the  moment  at  x  =  0  is  -Mo  and  the  moment  at 
x  =  0.58579*L  is  +Mo,  then  the  moment  can  not  exceed  Mo 
in  magnitude  anywhere  on  the  beam.  This,  then, 
constitutes  a  statically  admissible  moment 
distribution.  The  Lower  Bound  Theorem  essentially 
states  that  if  a  statically  admissible  moment 
distribution  can  be  found  that  satisfies  equilibrium 
and  the  boundary  conditions,  and  does  not  violate  the 
yield  criterion  (in  this  case,  M(x)<Mo  for  all  x),  then 
the  load  P  that  corresponds  to  this  moment  distribution 
is   a   lower   bound   to   the   exact   solution.   For  the 

L 

example   of   Figure   6,   the   value  of  P  =11.65685  is  a 

lower   bound,   but   in  this  case  it  also  corresponds  to 

u 
the   exact  solution.   If  the  upper  bound  load,  P  ,  were 

determined  at   the   same  point  on  the  beam,  the  same  P 

U   L 

would   result,  P  =P  ,  indicating  the  exact  solution  was 
obtained. 
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EFFECTIVE  PLASTIC  MOMENT  OF  A  STIFFENER 


Several  formulations  are  available  to  calculate  the 
effective  plastic  moment  of  a  given  I-beam,  yet  the 
problem  of  incorporating  a  general  formula  in  a 
computer  program  for  an  I-beam,  where  the  user 
specifies  the  thickness  and  breadth  of  the  upper  flange 
(plate  thickness  and  stiffener  spacing)  separately  from 
the  remaining  stiffener  scantlings,  presents  an 
interesting  subtlety.  The  program  must  be  able  to 
distinguish  between  the  case  where  the  neutral  exis 
lies  in  the  web  and  the  case  where  the  neutral  axis 
lies  in  the  plate.  A  simple,  accurate,  and  effective 
technique  to  accomplish  this  is  given  in  reference  {7}. 
Typical  stiffened  plate  geometry  is  shown  in  Figure  7. 
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Figure  7.    Stiffened  Plate  Geometry  for  Plastic  Section 
Modulus 


When   r  +r^>l,  the  plastic  neutral  axis  lies  in  the 
web  and  the  formula  for  the  plastic  section  modulus, Z, 

is: 

2 
Z  =  B  Dt  {r,  (r  +2r  +2)+2r  r^  (r  +l)-(r  -1)  }/4r,  . 
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When  r  +rp<l,  the  plastic  neutral  axis  lies  in  the 
upper  flange  (plate)  and  the  plastic  section  modulus  is 

given  by: 

2 
Z=B^Dt^{2r^+4r2+r^(l+2r^+2r2)-r^(r^+r2)  +2r2r4}/4. 

From  chapter  two,  the  value  of  effective  breadth  to 
be  incorporated  in  the  calculation  of  plastic  section 
modulus  was  shown  to  be  the  full  stiffener  spacing  {8}, 
without  an  arbitrary  cutoff  value.  CIRCPLAT  sets  B 
equal  to  stiffener  spacing  for  all  calculations. 


GRILLAGE  COLLAPSE  LOAD 


In  reference  {10},  Hodge  applied  the  Upper  Bound 
Theorem  to  a  rectangular  grid  (hereafter  called  a 
"grillage")  in  order  to  determine  its  collapse  load. 
In  the  method  development,  it  was  assumed  that 
concentrated  loads  were  applied  at  the  joints  (or 
nodes),  of  the  grillage,  although  the  method  could  be 
made  applicable  to  distributed  loads.  Having  assumed  a 
collapse   mode   and  determined  the  corresponding  upper 
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bound  collapse  load,  it  is  necessary  to  construct  a 
statically  admissible  bending  moment  distribution  to 
verify  that  the  assumed  collapse  mode  is  feasible.  If 
there  are  no  other  collapse  modes  that  would  yield  a 
lower  collapse  load  and  the  moment  distribution  is 
statically  admissible  (the  absolute  value  of  all 
moments  are  at,  or  lower  than,  the  plastic  yield 
moment),  then  the  correct  collapse  load  has  been 
determined.  Reference  {10}  contains  numerous  examples 
of  the  above  procedure. 

This  method  can  be  carried  out  with  little 
difficulty  for  regular  rectangular  grillages,  where  the 
loads  at  all  nodes  are  equal.  A  regular  rectangular 
grillage  is  one  where  the  spacing  is  uniform  between 
longitudinal  beams,  and  the  spacing  is  also  uniform 
between  transverse  beams.  In  cases  where  these 
requirements  are  not  met,  it  becomes  difficult  to 
predict  a  deformation  pattern  that  will  result  in  an 
improved  upper  bound.  Moreover,  the  assumed 
deformation  pattern  may  not  have  a  statically 
admissible  moment  distribution. 
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Reference  {5},  by  Hodge  and  Malone,  presented  a 
feasible  method  to  overcome  these  difficulties  in 
obtaining  a  correct  collapse  load.  The  proposed  method 
involves  taking  a  Lower  Bound  approach  using  Linear 
Programming.  The  load  associated  with  a  statically 
admissible  moment  distribution  is,  by  the  Lower  Bound 
Theorem,  a  lower  bound  to  the  true  collapse  load.  If 
the  load  can  be  maximized  while  the  moment  distribution 
is  maintained  statically  admissible,  the  maximum  lower 
bound  load  should  be  the  correct  collapse  load.  Linear 
Programming,  explained  in  Appendix  B,  provides  a 
simple,  computer  method  to  maximize  the  load,  subject 
to  the  constraints  of  static  admissibility  and  the 
equations  of  equilibrium. 

The  grillage  of  Figure  8  will  be  discussed  in 
detail  as  an  example  of  the  use  of  the  Lower  Bound 
Theorem  in  conjunction  with  Linear  Programming. 
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(a)  Grillage  with  load  (b)  Symmetric  node   (c)  Detail  of 

arrangement         Node  B 

Figure  8.    2x2  Grillage  with  Non-Uniform  Loading. 


For  this  problem  the  beams  are  simply  supported  at 
all  boundaries.  The  plastic  yield  moment  of  the 
horizontal  and  vertical  beams  is  Mo.  Moments  are 
positive  when  they  produce  compression  in  the  top 
fibers. 
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From  Figure  8b  we  can  see  that  there  is  symmetry 
about  a  line  through  nodes  A  and  C.  This  symmetry  can 
be  used  to  more  effectively  solve  the  problem.  As  a 
result  of  symmetry,  the  number  of  unknown  moments  is 
reduced  from  8  to  4.  The  numbers  on  the  ends  of  the 
beam  spans  indicate  the  unknown  moments  in  Figure  8b. 

Because  the  problem  of  plastic  collapse  is  only 
concerned  with  the  statically  admissible  moment  field, 
support  reactions  are  not  required  to  be  solved  for 
and,  as  a  result,  are  not  listed  as  unknowns. 

Since  the  loads  are  applied  at  the  nodes,  there  are 
no  midspan  loads  and  so  for  each  span,  the  shear  forces 
at  the  ends  must  be  equal  and  opposite.  Moment 
equilibrium  at  node  B  for  each  span  means  that 
V^=(M.-0)/L,  V2=(M^-M^)/L,  etc,  from  Figure  8c.  The 
sum  of  the  shear  forces  at  each  node  must  equal  the 
applied  load  at  that  node.   For  node  B,  the  result  is: 

(2M^+2M2-M^-M^  ) /L  =  0 
or  2M  +2M  -M, -M   =   0 

Applying  similar  reasoning  to  the  other  two 
independent  nodes: 
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4Mj^-2M2-P  =  0 
4M2j,-2Mo    =  0 


Notice   that   there   are   four  unknown  moments,  and 
only  three  independent  equilibrium  equations. 

In   addition   to   the  equations  of  equilibrium,  the 

linear   programming   problem   requires  the  equations  of 

constraint.    Each   of  the  moments  is  constrained  to  be 

equal   to   or   less   than   the   plastic   yield  moment, 

however,   this  moment  can  be  positive  or  negative.   The 

linear   programming   constraint   values   must   all   be 

positive.    In   order   to  meet  this  requirement,  a  new 

variable   is   created,  X. =M.  +Mo  i=l,2,3,4  so  that  while 

11 

-Mo<M.<+Mo,   0<X.<2Mo.   Therefore,  X.  is  constrained  to 

1        ~  1*  1 

positive   values.    The   collapse   load   is   the   fifth 
unknown,  called  X-. 

Substituting  M.  =X.  -Mo  into  the  three  equilibrium 

^11  ^ 

equations  the  result  is: 

-X^+2X2+2X  -X^  =  2Mo 

4Xi-2X      -X   =  2Mo 

-2X_+4X^  =  2Mo 
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The  load,  which  is  the  function  to  be  maximized  can  be 
defined  by  a  separate  equation  (called  the  objective 
function ) : 

X5=  P 

By  the  above  procedure,  all  of  the  requirements  for 
establishing  a  linear  programming  problem  have  been  met. 
The  next  step  is  to  format  the  information  to  facilitate 
input  to  the  Linear  Programming  routine.  The  required  input 
format  is  a  function  of  the  linear  program  to  be  utilized. 
At  Massachusetts  Institute  of  Technology,  one  of  the 
standard  linear  programs  available  in  the  International 
Mathematical  and  Statistical  Library  (IMSL)  is  called  ZX3LP 
(see  reference  {14}).  This  program  was  utilized  by  CIRCPLAT 
and  is  explained  fully  in  Appendix  B. 

ZX3LP  requires  the  input  to  be  arranged  in  3  matrices. 
Matrix  A  contains  the  coefficients  of  the  unknown  variables 
in  the  equations  of  constraint,  followed  by  the  coefficients 
of  the  unknown  variables  in  the  equilibrium  equations. 
Matrix  B  contains  the  values  from  the  right-hand  side  of  the 
constraint  equations,  followed  by  the  values  from  the  right 
hand  side  of  the  equilibrium  equations.  Matrix  C  contains 
the   coefficients   of  the  unknown  variables  in  the  objective 
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function.    Matrices   A,B,   and  C  for  the  sample  problem  of 
Figure  8  are  shown  in  Figure  9. 


A  Matrix 


B  Matrix 
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Figure  9. 


Linear  Program  ZX3LP  Input  Matrices. 
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The   primary   solution  of  the  ZX3LP  program  is  contained 
in   array   PSOL(N).    The   values   in   array  PSOL(N)  are  the 

values   of  X^,X2,Xo, ,Xj^.   For  the  problem  of  Figure  8, 

the  values  of  PSOL(l)  through  PS0L(4)  must  be  adjusted  by 
the  value  Mo  to  give  the  actual  nodal  moments.  This  is 
because  Mj^=Xj^-Mo.  The  value  of  the  collapse  load  is  given 
by  PS0L(5)  for  this  case. 


REPLACEMENT  THEOREM 


The  derivation  of  a  method  for  determining  the  collapse 
load  of  a  grillage  has,  up  to  now,  been  for  grillages  with 
concentrated  loads  at  beam  intersections.  In  extending  this 
method  to  encompass  the  problem  of  a  uniformly  loaded 
circular  plate,  it  is  necessary  to  discuss  the  correlation 
between  the  concentrated  load  problem  and  the  case  of  the 
uniform  load.  Reference  {10}  contains  a  discussion  and 
proof  of  the  Replacement  Theorem.  The  essence  of  the 
Replacement  Theorem  is  that  if  a  given  beam  span  of  constant 
plastic  yield  moment.  Mo,  is  acted  upon  by  a  distributed 
load  whose  direction  does  not  change  sign  on  that  span,  the 
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span  is  not  made  stronger  if  the  distributed  load  is 
replaced  by  equipollent  concentrated  loads.  That  is,  the 
collapse  load  of  a  span  under  a  set  of  concentrated  loads 
will  be  less  than  that  for  the  same  span  under  the 
equipollent  distributed  load. 

An  equipollent  force  system  is  defined  in  reference  {12} 
as  one  in  which: 

1.  The  resultant  of  the  forces  of  one  system  is  equal  to 
the  resultant  of  forces  of  the  other  system. 

2.  The  sum  of  the  moments  of  one  system  about  an  arbitrary 
base  point  is  equal  to  the  sum  of  the  moments  of  the 
other  system  about  the  same  point. 

Using  the  above  theorem,  if  the  uniform  load  on  the 
grillage  is  replaced  by  an  equipollent  system  of 
concentrated  forces  at  nodal  points,  then  the  method  of 
grillage  analysis  just  developed,  which  computes  a  "best" 
lower  bound  collapse  pressure  for  concentrated  forces,  can 
be  applied  directly  to  compute  a  conservative  "best"  lower 
bound  for  a  uniformly  distributed  load. 
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Looking  at  the  section  of  stiffened  plate  in  Figure  10, 
it  can  be  seen  that  the  pressure  P  distributed  over  an  area 
a  X  b  can  be  replaced  by  a  concentrated  force  of  magnitude 
P(a  X  b). 


Figure  10.    Section  of  a  Stiffened  Plate 


If  the  spacing  is  uniform  in  the  x-direction  and  is  also 
uniform  (not  necessarily  the  same)  in  the  y-direction,  then 
the  area  a  x  b  will  be  the  product  of  the  stiffener  spacing 
in  the  x  and  y  direction.  From  this,  it  can  be  seen  that 
for  a  regular  rectangular  grillage,  the  concentrated  force 
at  every  node  will  be  the  same.  In  the  case  of  a  grillage 
stiffened  circular   plate,  where  the  stiffeners  are  equally 
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spaced  along  the  main  diameters  (x  and  y  directions),  this 
same  replacement  process  works  exactly  the  same  for  nodes 
not  near  the  circumference  of  the  circle.  The  replacement 
process  for  nodes  near  the  circumference  is  not  exact,  as 
can  be  seen  in  Figure  11.  There  are  nodes  where  the 
replacement  force  will  be  too  high,  and  nodes  where  the 
replacement  force  will  be  too  low,  but  overall,  the  area  of 
the  circle  is  greater  than  the  area  associated  with  all  of 
the  nodes. 


Figure  11.    Top  View  of  a  Stiffened  Circular  Plate. 
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It  is  necessary  to  ensure  that  the  estimate  provided  by 
the  application  of  the  Replacement  Theorem  remains 
conservative  for  the  case  of  the  stiffened  circular  plate. 
In  CIRCPLAT,  the  forces  constrained  to  be  equal  at  all  of 
the  nodes.  Then,  after  the  collapse  load  is  determined,  the 
collapse  pressure  is  calculated  by  dividing  the  sum  of  the 
nodal  forces  by  the  total  area  of  the  circular  plate 
[P=(Sum  of  F)/A  .  Since  the  area  in  the  denominator  is 
greater  than  that  associated  with  the  nodal  forces  (the 
rectangle  of  Figure  11.),  the  pressure  will  be  conservative. 
The  accuracy  of  this  method  is  improved  as  the  number  of 
beams  in  the  grillage  is  increased. 


UNSUPPORTED  PLATE  COLLAPSE  LOAD. 


Calculating  the  collapse  load  of  the  grillage  does  not 
provide  any  information  about  the  possibility  of  collapse  of 
the  sections  of  plating  between  stiffeners.  To  be  complete, 
analysis  of  the  stiffened  plate  must  include  a  determination 
of  the  collapse  pressure  of  the  unsupported  plate  segments. 
Figure   12   indicates   the   variety   of   shapes  possible  for 
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unsupported  segments  of  a  stiffened  circular  plate.  In 
order  to  simplify  the  collapse  load  calculation,  use  is  made 
of  a  corollary  to  the  Lower  Bound  Theorem  discussed  in 
reference  {13}.  This  corollary  states  that  you  cannot 
decrease  the  collapse  load  of  a  structure  by  increasing  the 
strength  of  any  part. 


Figure  12.   Unsupported   Plate   Segments   of   Stiffened 
Circular  Plate. 


In  Figure  12  we  can  see  that  if  the  rectangular  sections 
of   unsupported   plating  are  a  x  b,  the  unsupported  sections 
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along  the  perimeter  of  the  circle  are  less  than  a  x  b. 
These  sections  can  be  considered  to  be  rectangles  of 
dimensions  a  x  b  with  a  wall  built  in  to  the  rectangle  for 
increased  support.  Here,  invoking  the  corollary  to  the 
Lower  Bound  Theorem  discussed  above,  it  can  be  seen  that 
these  rectangles  strengthened  by  walls  cannot  have  a  lower 
collapse  load  than  the  full  rectangular  section. 

Reference  {1}  contains  a  conservative  formula  for  the 
lower  bound  collapse  load  of  a  clamped  rectangular  plate 
subject  to  a  uniform  pressure.  The  unsupported  sections  of 
plate  are  best  modelled  as  clamped  because  the  rotation  at 
the  support  is  zero. 

This  formula,  implemented  in  CIRCPLAT,  is: 

P  =  (16Mo/b^)(l+S^) 
where  a=plate  length 
b=plate  breadth 
S=b/a 

This  equation  is  a  lower  bound  associated  with  the 
Johansen  yield  curve. 
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CHAPTER  4    CIRCPLAT  PROGRAM 


PROGRAM  DESCRIPTION  AND  USE 


CIRCPLAT  was  written  using  FORTRAN-77  programming 
language  on  a  VAX-11  computer  with  the  VAX/VMS  operating 
system.  The  program  was  designed  to  be  run  interactively 
with  a  minimum  of  required  information  from  the  program 
user,  while  providing  the  user  with  the  choice  of  inputting 
a  greater  amount  of  detailed  information  if  it  is  considered 
necessary. 

CIRCPLAT  was  designed  to  handle  a  maximum  of  twenty 
stiffeners  (horizontal  plus  vertical)  but  this  capacity 
could  be  increased  by  modifying  the  DIMENSION  statement  at 
the  beginning  of  the  program.  Similarly,  CIRCPLAT  can 
accommodate  up  to  a  maximum  of  12  stiffeners  in  any  one 
direction  (horizontal  or  vertical). 

The  stiffener  spacing  is  assumed  by  CIRCPLAT  to  be 
uniform  along  the  horizontal  and  vertical  diameters.   Thus, 
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when  the  user  inputs  the  total  number  of  vertical 
stiffeners,  the  total  number  of  horizontal  stiffeners,  and 
the  bulkhead  diameter,  the  program  computes  the  x  and  y 
coordinates  of  every  beam  intersection,  including  the 
coordinates  of  beam-boundary  intersections.  These 
coordinates  are  used  to  determine  the  span  lengths  that  are 
needed  for  calculation  of  the  local  shear  forces  in  the 
nodal  equilibrium  equations. 

CIRCPLAT  takes  advantage  of  the  symmetry  that  exists 
about  the  horizontal  and  vertical  diameters  and  reduces  the 
remaining  input  required  (and  the  calculations)  to  only 
those  concerned  with  the  upper  left  quadrant  of  the  circle 
(which  includes  any  beam(s)  that  lie  directly  on  the 
horizontal  and/or  vertical  diameters).  During  user  input, 
the  horizontal  beams  are  numbered  1,2, 3, etc.,  starting  from 
the  top,  counting  down  to  (and  including)  the  horizontal 
beam  at  the  centerline  (if  there  is  one).  Similarly, 
vertical  beams  are  numbered  1,2, 3, etc.  from  left  to  right. 

Inputs  for  the  material  yield  stress,  plate  thickness, 
then  the  stiffener  scantlings  are  requested  by  the  program 
and  the  plastic  yield  moments  are  calculated.  As  a 
convenience,  the  plastic  moment  resulting  from  the  first  set 

Chapter  4   CIRCPLAT  Program  49 


of  vertical  stiffener  scantlings  is  assigned  to  all  vertical 
stiffeners,  the  user  can  then  modify  this  value,  stiffener 
by  stiffener,  as  required.  The  user  has  the  option  of 
revising  the  plastic  moment  of  each  stiffener  by  either 

1.  Inputting  a  full  set  of  revised  scantlings  for  that 
stiffener  from  which  the  program  calculates  the  plastic 
yield  moment;  or 

2.  Inputting  the  plastic  yield  moment  value  directly  for 
the  specified  beam. 

Once  the  vertical  beams  have  been  fully  specified,  the 
identical  routine  is  run  for  the  horizontal  beams. 

The  option  of  specifying  the  plastic  yield  moment  value 
directly  for  an  individual  beam  was  added  to  allow  the  user 
to  insert  a  calculated  plastic  moment  corrected  for  the 
effect  of  shear  stress  or  axial  force  if  it  should  become 
necessary.  The  method  of  obtaining  these  corrections  is 
discussed  in  detail  in  reference  {4}. 

CIRCPLAT  was  designed  for  analysis  of  a  grillage  where 
all   beams   are   prismatic.   For  those  cases  involving  beams 
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that  vary  in  cross-section  along  their  length,  the  minimum 
plastic  yield  moment  should  be  used  to  obtain  a  conservative 
failure  load,  or,  the  program  itself  will  have  to  be 
modified  to  allow  for  changes  in  stiffener  geometry  along  a 
given  beam  length. 

For  those  beams  where  shear  deformations  or  effects  must 
be  taken  into  account,  if  the  corrected  plastic  yield  moment 
is  less  than  the  minimum  plastic  yield  moment  for  the 
uncorrected  beam,  the  corrected  value  should  be  input. 
Otherwise,  to  remain  conservative,  the  minimum  uncorrected 
plastic  yield  moment  should  continue  to  be  used. 

The  current  version  of  CIRCPLAT  is  designed  to  analyze 
the  built-in  submersible  bulkhead.  The  program  treats  the 
grillage  as  being  fully  clamped  at  all  boundaries. 

At  this  point,  the  user  input  is  completed,  and  the  program 
begins  to  calculate  the  parameters  that  it  will  have  to 
input  to  the  ZX3LP  linear  programming  subroutine. 
Basically,  the  remainder  of  the  program  is  designed  to 
accomplish  four  functions: 
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1.  Loading  the  correct  values  into  the  proper  locations  in 
Matrices  A,B,and  C  as  indicated  in  the  example  of 
Chapter  Three. 

2.  Establishing  the  exact  matrix  DIMENSION  values  required 
by  ZX3LP. 

3.  Calculating  the  failure  load  of  the  unsupported  sections 
of  plating  between  stiffeners. 

4.  Formatting  the  solution  produced  by  ZX3LP  to  produce  a 
properly  labelled  output. 

The  solution  is  stored  automatically  in  a  data  file 
labelled  FOR018.DAT.  This  data  file  contains  a  listing  of 
the  input  data  and  the  output  data. The  output  data  consists 
of  a  listing  of  the  following  items: 

1.  The  failure  load  of  the  weakest  segment  of  unsupported 
plating. 

2.  The  failure  load  of  the  grillage  as  a  whole.  This  value 
is  unrelated  to  the  value  in  (1)  above  (i.e.  they  are 
treated  as  separate  problems). 
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3.  The  normalized  value  (with  the  corresponding  plastic 
moment  value  used  to  normalize  it)  of  the  moment  present 
at  failure  at  each  of  the  vertical  beam  nodes.  This  is 
followed  by  a  similar  listing  for  all  horizontal  nodes. 
The  listing  of  node  moments  is  only  for  those  nodes  in 
the  upper  left  quadrant  of  the  plate.  The  rest  are 
symmetric.  Also,  the  vertical  nodes  are  numbered  from 
the  top  of  the  beam  downward,  working  from  left  to 
right.  In  a  similar  fashion,  the  horizontal  nodes  are 
listed  from  left  to  right,  working  from  top  to  bottom. 
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CHAPTER  5    RESULTS  AND  DISCUSSION 


COMPARISON  OF  RESULTS 


The  results  of  CIRCPLAT  were  compared  to  the  results 
obtained  from  an  actual  scale  model  failure  test  of  a 
stiffened  circular  bulkhead  (Reference  {16}).  The  ratio  of 
CIRCPLAT  predicted  collapse  pressure  to  experimental 
collapse  pressure  was  0.61.  There  are  a  number  of  valid 
reasons  why  the  CIRCPLAT  prediction  was  39%  lower  than  the 
experimental  value. 

First,  the  horizontal  and  vertical  stiffeners  on  the 
bulkhead  under  test  were  not  prismatic.  That  is,  the 
scantlings  varied  significantly  along  the  stiffener  length. 
To  accommodate  possibly  large  shear  effects  on  the  main 
horizontal  girder,  the  web  depth  increased  in  a  linear 
fashion  (over  approximately  one  quater  span)  at  each  end  to 
a  web  depth  almost  twice  that  of  the  center  span. 
Variations  of  web  depth  of  a  similiar  magnitude  were  also 
common   on  vertical  stiffeners.   Because  the  program  assumes 
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prismatic   beams,   and   to   ensure  conservative  results,  the 
smallest  value  for  web  depth  was  input  for  every  stiffener. 

Secondly,  the  bulkhead  design  tested  incorporated  both 
HY-80  and  HY-130  steels.  Again,  because  the  program  allowed 
for  only  one  material  yield  strength,  and  to  ensure 
conservative  results  only  the  lower  yield  stress  was  used  in 
the  analysis. 

Finally,  a  less  important  factor,  the  program  can 
presently  handle  a  maximum  of  twelve  stiffeners  in  any  one 
direction.  There  are  twelve  tertiary  horizontal  stiffeners 
and  one  main  horizontal  girder  for  the  tested  bulkhead  (for 
a  total  of  thirteen).  To  maintain  the  odd-symmetry 
relationship,  and  keep  the  total  under  twelve,  it  was 
necessary  to  reduce  the  number  of  tertiary  stiffeners  to 
ten.  This  has  a  minor  effect,  as  indicated  by  the  fact  that 
the  effect  of  tertiary  stiffeners  was  ignored  entirely  in 
the  analysis  of  reference  {4}.  Again,  the  results  will  be 
more  conservative. 

As  a  result  of  the  above  three  factors,  CIRCPLAT  results 
can  be  assumed  to  be  very  conservative;  this  is  borne  out  by 
the  above  comparison. 
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The  previous  results  were  for  a  strengthened  version  of 
an  existing  circular  bulkhead  design.  For  the  original  (not 
strengthened)  version  of  this  bulkhead,  scantling  data  was 
readily  available,  but  the  given  collapse  load  data  was 
based  on  design  values,  rather  than  actual  failure  testing. 
CIRCPLAT  was  run  and  results  were  compared  to  the  design 
collapse  load  for  the  circular  bulkhead  of  the  original 
design.  When  CIRCPLAT  was  run  with  tertiary  stiffeners  not 
included,  the  ratio  of  CIRCPLAT  collapse  pressure  to  design 
collapse  pressure  was  0,83.  When  CIRCPLAT  was  run  with  all 
tertiary  stiffeners  (8  of  them)  included,  the  ratio  was 
1.27.  The  correlation  here  is  considered  good,  and  the 
effect  of  tertiary  stiffeners  for  this  case  would  seem  to  be 
significant.  Failure  load  analysis  of  the  unsupported 
plating  segments  for  the  original  design  indicated  that  the 
unsupported  plating  was  the  weakest  portion.  The  ratio  of 
theoretical  collapse  pressure  to  design  collapse  pressure 
for  the  unsupported  plate  (clamped  rectangular  plate)  was 
0.73  for  an  Upper  Bound  solution,  and  0.49  for  a  Lower  Bound 
solution.  This  indicated  the  bulkhead  plating  would  not 
withstand  full  design  pressure.  This  problem  did  not  show 
up  in  the  improved  bulkhead  design. 
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In  reference  {10},  Hodge  provides  a  formula  for  the 
exact  collapse  load  of  a  uniform  circular  plate,  with  a 
fully  clamped  boundary.   The  formula  reduces  to 

P  =  2.8140^  t^/R^ 
where  P  =  uniform  collapse  pressure 
(To   =  material  yield  stress 
t  =  plate  thickness, 
R  =  Radius  of  plate. 

For  a  plate  10  feet  in  diameter,  2  inches  thick,  with  a 
yield  stress  of  50,000  psi,  the  exact  uniform  collapse 
pressure  is  156  psi. 

Referring  to  the  formulas  of  Chapter  3  for  calculating 
the  plastic  section  modulus  of  a  stiffened  plate,  CIRCPLAT 
can  be  used  to  accurately  model  a  circular  plate  as  a 
grillage  by  setting  all  web  and  flange  scantlings  equal  to 
zero,  except  for  the  web  depth.  The  web  depth  can  be  any 
arbitrary  value,  as  it  factors  out  of  the  plastic  section 
modulus   formula,   but  it  must  not  be  zero,  because  it  is  in 

the  denominator  of  r   and  r,  and  would  cause  these  values  to 

J  '4- 

become  singular.  For  convenience,  input  web  depth  to 
CIRCPLAT  as  1.0,  and  input  all  other  web  and  flange 
scantlings  as  zero.   As  a  result,  the  plate  will  be  modelled 
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by  CIRCPLAT  as  a  grillage  of  stiffeners  with  rectangular 
cross  section  (bxh)  where  b  is  stiffener  spacing,  determined 
by  the  number  of  stiffeners  the  user  inputs  (more  stiffeners 
will  give  better  results);  and  h  is  plating  thickness. 

For  the  same  plate  analyzed  above  using  Hodges' s 
formula,  CIRCPLAT  calculated  a  failure  load  of  154.7  psi, 
only  1%  less  than  the  exact  value.  The  plate  was  modelled 
for  this  case  as  ten  vertical  and  ten  horizontal  stiffeners. 
This  result  is  significant  in  that  it  points  out  that 
CIRCPLAT  does  yield  a  conservative  yet  very  close,  lower 
bound  solution,  in  a  case  where  the  input  values  can  closely 
approximate  the  exact  physical  problem,  and  also  verifies 
that  a  uniform  load  can  be  approximated  quite  well  by  point 
loads  at  all  nodes. 


SUGGESTED  IMPROVEMENTS  TO  CIRCPLAT 


In  order  to  reduce  the  amount  of  conservatism  indicated 
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by  the  comparison  of  CIRCPLAT  results  to  actual  stiffened 
plate  failure  test  data,  several  improvements  to  the  program 
are  recommended  for  implementation: 

1.  Modify  the  input  routine  to  provide  the  user  with  the 
ability  to  input  a  different  set  of  horizontal  and/or 
vertical  stiffener  scantlings  at  each  nodal  point.  This 
would  enable  CIRCPLAT  to  more  closely  approximate  a 
non-prismatic  beam; 

2.  Modify  the  input  routine  to  provide  the  user  with  the 
ability  to  input  a  different  yield  stress  value  for  each 
stiffener,  or  for  the  circular  plate.  This  would, 
again,  provide  the  user  with  more  flexibility  in 
modelling  a  real,  physical  problem; 

3.  Modify  the  program  to  include  a  variable  "fixity" 
factor,  f,  at  boundary  nodes,  ranging  in  value  from  zero 
to  one.  This  would  allow  the  user  to  model  a  simple 
support  (f=0.),  a  fixed  support  (f=l.),  or  some  value  in 
between . 
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CONCLUSION 


In  conclusion,  it  is  felt  that  the  linear  programming 
method  of  grillage  analysis,  based  on  the  Lower  Bound 
Theorem,  provides  a  simple,  yet  powerful  tool  for  accurately 
determining  the  collapse  load  of  a  given  structure. 
CIRCPLAT  was  designed  for  the  problem  of  the  circular  plate 
stiffened  by  a  grillage  network,  with  fixed  (clamped) 
boundaries.  However,  the  applicability  of  this  technique 
can  be  extended  very  easily  to  all  types  of  grillage 
collapse  analysis.  It  is  strongly  felt  that  this  analysis 
tool  yields  results  that  are  consistent  with  empirical  data. 
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APPENDIX  A   PROGRAM  LISTING 


Start 
IE 


INPUT 

No.  of  Vert.  Stiffeners 
No.  of  Horz.  stiffeners 
Plate  Diameter 
Plate  Thickness 
Material  Yield  Stress 


INPUT 

Scantlings  of  most  common 
Vertical  Stiffener 
(Web  Thickness/Depth, 
Flange  Thickness/width ) 


YES 


NO 


Repeat  process  for 
Horizontal  Stiffeners 


INPUT  1 ,or  2, 
above 


-® 
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Calculate  Failure 
Load  of  Unsupported 
Plate 


call  Subroutine  ZX3LP 
to  solve  Linear  Programming 
Problem 


Store  OUTPUT  in  file 
F0R018.DAT 


END 
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10  C23456789 

20  DIMENSION  XCOORD ( 8 , 8) , YCOORD ( 8 , 8 ) , VERTMO ( 8 ) , 

30  1  H0RZM0(8)  ,MV(8,8) ,MH(8, 8)  ,NNODE(40)  ,F(40)  , 

40  1  NVERT(40)  ,MLIN( 80)  ,MMAX(80)  ,C( 80)  ,AA( 120 , 80)  , 

50  1  B(  100)  ,HORZL(40)  , VERTL(40)  ,KHORZ(40)  ,KVERT(40)  , 

60  1  MPV(45)  ,MPH(45)  ,NVLEN(45)  , IFLAG(45)  ,PSOL( 102 ) , 

70  1  A(120,80) ,DSOL( 120) ,RW(8000) ,IW(304) 

80 

90  REAL  HORZMO,MV ,MH,MMAX,MLIN,HORZL,MPV ,MPH,HFLAG,MP 

100 

110  c 

120  C    USER  INTERACTIVELY  INPUTS  BEAM  AND  PLATE  PARAMETERS  FOR 

130  C    THE  STIFFENED  CIRCULAR  PLATE  OF  INTEREST.  SYMMETRY  ABOUT 

140  C    THE  HORIZONTAL  AND  VERTICAL  CENTERLINES  IS  ASSUMED. 

150  C    INPUT  PARAMETERS  AND  CALCULATIONS,  THEN,  ARE  FOR  THE 

160  C    UPPER  LEFT  QUADRANT  OF  THE  CIRCULAR  PLATE. 

170  C 
180 

190  C234567 

200  PRINT  *,  '   THIS  PROGRAM  WAS  DESIGNED  FOR  A  LIMIT  OF' 

210  PRINT  *,'   TWENTY  TOTAL  STIFFENERS  ( HORI Z . +VERT . ) ' 

220  PRINT  *,'   TWELVE  MAX.  IN  ONE  DIRECTION.' 

230  Y=SKIP(5) 

240  PRINT  *,  '   THIS  PROGRAM  ASSUMES  THAT  HORIZONTAL  AND' 

250  PRINT  *,  '   VERTICAL  STIFFENERS  ARE  EQUALLY  SPACED  ALONG' 

260  PRINT  *,  '   THE  HORIZONTAL  AND  VERTICAL  CENTERLINES,  AND' 

270  PRINT  *,  '   ARE  T-STIFFENERS  WITH  AN  EFFECTIVE  PLATE  WIDTH' 

2  80  PRINT  *,  '   MAKING  THEM  I -ST IFFENERS . ' 
290           Y=SKIP(5) 

300  PRINT  *,  '   ENTER  THE  NO.  OF  VERTICAL  STIFFENERS' 

3  10  READ  *,  N2 
320  Y=SKIP(5) 

330  PRINT  *,  '   ENTER  THE  NO.  OF  HORIZONTAL  STIFFENERS' 

3  40  READ  *,  M2 

350  Y=SKIP(5) 

360  PRINT  *,  '   ENTER  THE  DIAMETER  (FEET)  OF  THE  CIRC.  PLATE' 

370  READ  *,  DIAM 

380  D=DIAM*12. 

390  R=D/2. 

400  WRITE( 18,4000) 

410  4000  FORMATC IX, ' INPUT  INFORMATION'/) 

420  C234567 

430  WRITE( 18,4001 )  N2 

440  4001  F0RMAT( IX, 'NO.  OF  VERTICAL  ST IFFENERS :', 3 X , 13  /) 

450  WRITE(18, 4002 )  M2 

460  4002  FORMAT( IX, 'NO.  OF  HORIZONTAL  ST  I FFENERS :  '  ,3 X  ,  13  /) 

470  WRITE( 18, 4003 )  DIAM 


480     4003  FORMATC IX, 'PLATE  DI AMET ER :  '  , 3 X  ,  F6  .  2  ,  ' ( FEET  )  '  /) 

490  C 

500  C   CALCULATE  THE  X  (HORIZ.)  AND  Y  (VERT.)  COORDINATES  OF 

510  C   NODES  (BEAM  INTERSECTIONS  OR  BOUNDARIES) 

520  C 

53  0  C 

540  YDIST=D/REAL(M2+1 ) 

550  XDIST=D/REAL(N2+1) 

560  C  CHECK  ODD  OR  EVEN  SYMMETRY 

570  VAR1=REAL(N2 )/2.-IFIX(REAL(N2)/2.)-.l 

5  80  VAR2=REAL(M2 ) /2  . -IFIX( REAL(M2 ) H  .)-  .1 

590  IODDV=0 

600  IODDH=0 

610  IF( VARl .GT.O)  I0DDV=1 

620  IF(VAR2 .GT.O)  I0DDH=1 

630  IF(IODDV.EQ.O)  GO  TO  50 

640  KCOL=(N2+l)/2+l 

650  VFLAG=0.0 

660  GO  TO  51 

670     50    KC0L=N2/2+l 

680  VFLAG=1.0 

690     51    IF( lODDH.EQ.O)  GO  TO  55 

700  KROW=(M2+l)/2+l 

710  HFLAG=0.0 

720  GO    TO    56 

730  55  KROW=M2/2+l 

740  HFLAG=1.0 

750     56    CONTINUE 

760  C 

770  C   KCOL  =  NO.  OF  NODES  ON  LONGEST  HORIZ.  BEAM 

780  C   KROW  =  NO.  OF  NODES  ON  LONGEST  VERT.  BEAM 

790  C   (THIS  IS  INCLUDING  BOUNDARY  NODES  AND  IS  FOR 

800  C   THE  UPPER  LEFT  QUADRANT  ONLY) 

810  C   H/VFLAG  =  0  IF  ODD  SYMMETRY 

820  C  =  1  IF  EVEN  SYMMETRY 

83  0  C 

840  C234 

850  C   INITIALIZE  ALL  COORD.  MATRICES  TO  ZERO 

860  C 

870  DO  70  1=1, KROW 

880  DO  60  J=l ,KCOL 

890  XCOORD(I  ,  J)=0.0 

900  YCOORDd,  J)=0.0 

910     60        CONTINUE 

920     70    CONTINUE 

930  C 

940  C   FILL  IN  COORD.  MATRICES 

950  C 
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960  DO  110  NR0W=1,KR0W 

970  DO  100  NCOL=KCOL, 1 ,-1 

980  IF(NROW.GT.  1)  GO  TO  115 

9  90  XVAL=R-(REAL(KC0L-NC0L)+VFLAG/2 . )*XDIST 

1000  YVAL=SQRT(R*R-(R-XVAL)**2 ) 

1010  YMIN=(REAL(KROW)-2+HFLAG/2 . )*YDIST 

1020  TEST1=YVAL-YMIN 

1030  IF(TEST1  .LT.O)  GO  TO  110 

1040  XCOORD(NROW,NCOL)=XVAL 

1050  YCOORD(NROW,NCOL)=YVAL 

1060  GO  TO  100 

^070  115            XMAX=SQRT(R*R-((REAL(KR0W-NR0W)+HFLAG/2 .)*YDIST)**2 ) 

10  80  YMAX=SQRT(R*R-( ( REAL ( KCOL-NCOL) +VFLAG/2 . ) *XDI ST ) **2 ) 

1090  C234567 

1100  IF(NCOL.EQ. 1)  GO  TO  116 

1110  XVAL=R-(REAL(KC0L-NC0L)+VFLAG/2 . )*XDIST 

1120  XCOORD(NROW,NCOL)=XVAL 

1130  TEST2=XMAX-(R-XVAL) 

11^0  IF(TEST2  .LE.O)  GO  TO  105 

1150  106            YVAL=(REAL(KR0W-NR0W)+HFLAG/2  .  )*YDIST 

1160  TEST3=YMAX-YVAL 

1170  IF(TEST3 .LT.-YDIST)  GO  TO  110 

1180  IF(TEST3  .LE.O)  GO  TO  108 

1190  YCOORD(NROW,NCOL)=YVAL 

1200  GO  TO  100 

1210  105           XCOORD(NROW,NCOL)=R-XMAX 

1220  GO  TO  106 

1230  108           YCOORD(NROW,NCOL)=YMAX 

1240  GO  TO  110 

1250  116           XCOORD(NROW,NCOL)=R-XMAX 

^2^0  YC00RD(NR0W,NC0L)  =  (REAL(KR0W-NR0W)+HFLAG/2  .)*YDIST 

1270  GO  TO  110 

1280  100      CONTINUE 

1290  110   CONTINUE 

1300  C 

13  10  C   USER  INPUTS  STIFFENER  SCANTLINGS  TO  DETERMINE  PLASTIC  MOMENTS 

1^30  PRINT  *,  '    INPUT  MATERIAL  YIELD  STRESS  (PSI)  AND  PLATE' 

1340  PRINT  *,  '   THICKNESS  (INCHES)  AS  FOLLOWS:' 

1350  PRINT  *,  '   SIGMA-Y,  THICKNESS' 

1360  READ  *,  SIGMAY,T2                                                   ' 

1370  Y=SKIP(5) 

1380  PRINT  *,  '   THE  INITIAL  SCANTLINGS  YOU  INPUT  WILL  BE  ASSIGNED' 

}?^°  ^^INT  *,  '   TO  ALL  VERTICAL  STIFFENERS.  YOU  WILL  THEN  HAVE  THE 

J;°°  PRINT  *,  '   OPTION  OF  REVISING  SPECIFIC  STIFFENERS  BY  NUMBER  ' 

IflO  PRINT  *,  '   (BY  ENTERING  EITHER  SCANTLINGS  OR  PLASTIC 


1420  Y=SKIP(5) 

1430  IVERT=0 


MOMENT) 
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^^^°  PRINT  *,  '    INPUT  INITIAL  VERTICAL  BEAM  SCANTLINGS  AS  FOLLOWS 

1450  15    PRINT  *,  '  ENTER  WEB  THICKNESS , WEB  DEPTH,  BOTTOM  FLANGE  ' 

1^^°  PRINT  *,  '  THICKNESS, BOTTOM  FLANGE  WIDTH  (ALL  IN  INCHES)' 

1^70  READ  *,  T3  ,DW,T1  ,B1 

1^80  IF(IVERT.EQ.O)  GO  TO  4005 

1490  WRITE(18,4006)  NBE AM- 1 , T3  , DW , T 1  ,  B 1 

1500  4006  FORMATdX, 'VERTICAL  BEAM  NO  .  :  '  ,  3  X  ,  13  ,  / 

1510  11X,'WEB  THICKNESS: ' ,3X,F8.3 , '  INCHES'  / 

1520  11X,'WEB    DEPTH   :',3X,F8.3,'  INCHES',/ 

1530  IIX, 'BOTTOM  FLANGE  THICKNESS :', 3 X , F8 . 3 , '  INCHES'  / 

1540  IIX, 'BOTTOM  FLANGE    WIDTH   :',3X,F8.3,'  INCHES''//) 

1550  GO  TO  16 

1560  4005  WRITE(18,4004)  S IGMAY , T2  ,  T3  ,  DW , T  1  ,  B  1 

1570  4004  FORMATdX,'  VERTICAL  BEAMS  ',/ 

1580  IIX, 'MATERIAL  YIELD  STRES S : ' , 3 X , F 12 . 2 , '  PSI'  / 

1590  IIX, 'PLATE  TH ICKNES S : ' , 3 X , F 8 . 3 , '  INCHES'  / 

1^00  IIX,  'COMMON  WEB  THI CKNES S :  '  ,3 X , F 8 .3  ,  '  INCHES'  / 

1610  IIX, 'COMMON  WEB  DEPTH     :',3X,F8.3,'  INCHES''/ 

1^20  IIX, 'COMMON  BOTTOM  FLANGE  THI  CKNES  S  :  '  ,  3  X  ,  F  8  .  3  ',  '  INCHES',/ 

l^^Q  UX,  'COMMON  BOTTOM  FLANGE  WIDTH      :'  3X  F8  3''  INfHFc;''//^ 

1640  16  IF(IVERT.EQ.l)  GO  TO  21                ,-iX,t».J,   INCHES  ,//) 

1650  DO  20  1=2 ,KCOL 

1660  VERTM0(I)=MP(DW,XDIST,T2 ,B1,T1,SIGMAY,T3) 

1670  20    CONTINUE 

1680  GO  TO  22 

1^90  21    VERTM0(NBEAM)=MP(DW,XDIST,T2 ,B1 ,T1 ,SIGMAY  T3 ) 

1700  22     PRINT  *,  '   ARE  THERE  ANY  VERT.  BEAMS  WHICH  ARE  DIFFERENT' 

I'^IO  PRINT  *,  '   FROM  THOSE  INPUT  ABOVE  ■?  ' 

1720  PRINT  *,  '   0=N0,  1=YES' 

1730  READ  *,  IVERT 

1740  IF(IVERT.EQ.O)  GO  TO  30 

1750  Y=SKIP(2) 

1^50  PRINT  *,  '   SINCE  SYMMETRY  IS  ASSUMED,  INPUT  ONLY  THOSE  VERT  ' 

;;;°  P^I^NT  *,  -   beams  on  or  to  the  left  of  VERT.  CENTERLINE.' 

1780  Y=SKIP(2) 

^^^^  PRINT  *,  '   COUNTING  IN  SEQUENCE  LEFT  TO  RIGHT   INPUT  THE' 

Jf?^  PRl^NT  *,  '   BEAM  NO.  FOR  WHICH  THERE  ARE  NEW  SCANTLINGS' 

1810  READ  *,  NBEAM 

1820  NBEAM=NBEAM+1 

f^^^  PRl^NT  *,  '   DO  YOU  WANT  TO  INPUT  SCANTLINGS,  OR  PLASTIC  MOMENT 

840  PRINT  *,  '   0=INPUT  PLASTIC  MOMENT , 1 = INPUT  SCANTLINGS' 

1850  READ  *,  NOPT 

I860  IF(NOPT.EQ.  1)  GO  TO  15 

l^^O  PRINT*,  '    INPUT  PLASTIC  MOMENT' 

1880  READ  *,  VERTMO(NBEAM) 

1890  Y=SKIP(2) 

1900  GO  TO  22 

1910  30   IHORZ=0 
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920  PRINT  *,  '   NOW  A  SIMILAR  ROUTINE  WILL  BE  FOLLOWED  FOR' 

\ll^  PRINT  *,  '   HORIZONTAL  STIFFENERS.' 

950  .S    ll.ll    !'  '   '""'"^  "°^''*  "'^^  SCANTLINGS  AS  FOLLOWS:' 

960  PR  NT  *'   '   THTr^T^L'  ^"IC'^NESS,  WEB  DEPTH,  BOTTOM  FLANGE' 

1980  ReId  */t3,Dw!t?'bi    '  ''^^  ''^'''  "'"""  ^^^^  ^^  ^^^»^S)' 

^990  IF(IHORZ.EQ.O)  GO  TO  4007 

2000  WRITE(18,4008)  MBEAM-1 , T3 , DW , T 1  Bl 

2010  4008  FORMATdX, 'HORIZONTAL  BEAM  NO.-'  3X  13  / 

2020  11X.;WEB  THICKNESS:', 3X,F8. 3,'  INCHES'./ 

2040  v'^^^    '^^P^"   :',3X,F8.3,'  INCHES'  / 

2^50  v'-^nJ^''  "'^^''^'^  THICKNESS:', 3X,F8. 3  '  INCHES',/ 

2060  'go'to  34     '"''"    ''""   :  '  ,  3  X  ,  F  8  . 3  ,  '  INCHES  '  ',  )  /  ) 

2070  4007  WRITE( 18,4009)  T3 , DW  Tl  Bl 

2080  4009  FORMATdX,  'HORIZONTAL  BEAMS'  / 

2?00  IfX, 'COMMON  WEB  THICKNESS:', 3i,F8. 3,'  INCHES',/   ' 

2nO  X,  COMMON  WEB    DEPTH   :  '  ,  3  X ,  F  8 .3  ,  '  INCHES '  / 

2110  X,  COMMON  BOTT.FLG. THICKNESS:', 3X,F8. 3.'  INCHES',/ 

;J?n  o,^^^'  COMMON  BOTT.FLG.   WIDTH   :'  3X  F8  3  '  INCHF^'  //^ 

2130  34  IF(IHORZ.EQ.l)  GO  TO  31          ,  J  a  ,  j.  o .  j  ,   INCHES  ,//) 

21^0  DO  33  1=2 ,KROW 

2J6O  33    ,,/,',llf^^^--^'^^^^^^^^^^^2,Bl,Ti,SlGnAY,,3) 

2170  GO  TO  32 

lllo  I?           ^°?^!!°i^^?^^)=^P(DW,YDIST,T2,Bl,Tl,SIGMAY,T3) 

2200  ll.ll    !'     ""^^    ™^^    ^^^  ^°^I^-  «E^^S  WHicH  ARE  ' 

,,,"  PRINT  *,      DIFFERENT  FROM  THOSE  INPUT  ABOVE  ^' 

2220  ^^^^^  *'  '   °^^°'  1=YES' 

2230  READ  *,  IHORZ 

2240  IF(IHORZ.EQ.O)  GO  TO  40 

2250  Y=SKIP(2) 

2260  C234567 

2280  imi    !'  ^'   ^'"""^^  SYMMETRY  ASSUMED.  INPUT  ONLY  THOSE' 

2290  Y  =  SKIP(2)    "°^^^*  ^^^"^^    °^  °^  ^^°^^  ^°^^^-  CENTERLINE' 

23?0  PPTMT  *'  '   ^^"NTING  IN  SEQUENCE  FROM  TOP  TO  BOTTOM,' 

2320  ReId  */mBEAm''"'  ^  ''"'  '^  *  '°'  ^'^  ^^^  SCANTLINGS.' 

2330  MBEAM=MBEAM+1 

2350  PRINT  *'  '   n^TMPn/f  "'  '"'""  SCANTLINGS  OR  PLASTIC  MOMENT? 

2360  ReId  */mOPT                '  ''°^'^''  '  =  '^'""  SCANTLINGS' 

2370  IF(MOPT.EQ.l)  GO  TO  35 

23  80  PRINT  *,  '   INPUT  PLASTIC  MOMENT' 

^3  90  READ  *.  HORZMO(MBEAM) 

2400  Y=SKIP(2) 
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2410  GO  TO  32 

2420     40    CONTINUE 

2430    C 

2440    C 

2450    C   INITIALIZE  ALL  MV  AND  MH  FIELDS  TO  ZERO 

2460    C 

2470    C 

2480  DO  90  I=1,KR0W 

2490  DO  80  J=1,KC0L 

2500  MV(I,J)=0.0 

2510  MH(I,J)=0.0 

2520     80         CONTINUE 

2530     90    CONTINUE 

2540    C234567 

2550    C 

2560    C 

2570    C   INITIALIZE  NODAL  EQUATION  COUNTER  TO  ZERO  AND  COUNT 

2580    C   INTERNAL  NODES  FOR  EACH  LINE.  SUM  INTERNAL  NODES  TO 

2590    C   COUNT  NODAL  EQUATIONS 

2600    C 

2610    C 

2620  NODEEQ=0 

2630  DO  120  J=2 ,KR0W 

2  6  40  XMAX  =  SQRT(R*R-( ( RE AL ( KROW- J )+HFLAG/2  .  ) * YD  1ST  )  * *2  ) 

2  6  50  NNODE( J)=IFIX(XMAX/XDIST-VFLAG/2 .+0. 9  9  99) 

2660  NODEEQ=NODEEQ+NNODE(J) 

2670     120   CONTINUE 

2680    C 

2690    C 

2700    C   CREATE  A  ONE  DIMENSIONAL  ARRAY  OF  UNIT  NODAL  GRID  LOADS 

2710    C 

2720    C 

2  730  DO  130  K=l  .NODEEQ 

2740  F(K)=1 

2750     130   CONTINUE 

2760  KKK=1 

2770  IF(KKK.GT.O)  IBC=2 

2780  IF(KKK.GT.O)  GO  TO  13  1 

2790  PRINT  *,  '   ENTER  BOUNDARY  CONDITIONS:' 

2800  PRINT  *,  '    1=SIMPLE  SUPPORT,  2=FIXED  SUPPORT' 

2810  READ  *,  IBC 

2820  131       CONTINUE 

2  83  0  C 

2840         C 

2850    C    END  OF  USER  INPUT 

2860    C 

2870  WRITE(18,4014) 

2880     4014  F0RMAT(1X,'  OUTPUT  INFORMATION  ') 
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2890  IF(NOPT.EQ. 1 .AND. MOPT.EQ. 1 .AND. T3  .LT.O  .00001  )  GO  TO  4020 

2900  BB=YDIST/2. 

2910  L=XDIST/2. 

2920  IF(XDIST.GE. YDIST)  BB=XDIST/2. 

2930  IF(XDIST.GE. YDIST)  L=YDIST/2. 

2  940  PL  =  SIGMAY*(T2**2 / BB**2 ) * ( 1 +BB**2 /L**2  ) 
2950           WRITE( 18,4015)  PL 

2960  4015  F0RMAT(1X,'  UNSUPP  .  PLATE  SECTION  COLLAPSE  LO AD :  '  , F  1  0 . 4 ) 

2970  C 

2980  C 

2990  C   ASSIGN  MOMENT  FLAGS  TO  HORIZ.  BOUNDARY  NODES  lAW  INPUT  B.C. 

3000  C 

3010  C 

3020  C 

3030  C234567 

3040  4020  DO  140  NN=2 , KROW 

3050  MC0L=KC0L-NN0DE(NN) 

3060  MH(NN,MC0L)=1 . 

3070  IF(IBC.EQ.l)  MH ( NN , MCOL ) =0 . 0 

3080  140   CONTINUE 

3090  C 

3  100  C   COUNT  INTERNAL  NODES  ON  VERTICAL  BEAMS 
3110  C 

3 120  DO  150  L=2 ,KC0L 

3 13  0  YMAX=SQRT(R*R-( ( RE AL ( KCOL-L )+VFLAG/2 . ) *XDI ST ) **2 ) 

3 140  NVERT(L)=IFIX(YMAX/YDIST-HFLAG/2 .+0.9999) 

3150  150   CONTINUE 

3160  C234567 

3  170  C 

3  180  C   ASSIGN  MOMENT  FLAGS  TO  VERT.  BOUNDARY  NODES  lAW  INPUT  B.C. 

3190  C 

3200  C 

32 10  DO  160  JK=2 ,KC0L 

3220  MR0W=KR0W-NVERT( JK) 

3230  MV(MR0W, JK)=1 . 

3240  IF(IBC.EQ.l)  MV ( MROW , JK ) =0 . 0 

3250  160   CONTINUE 

3260  C 

3270  C   SET  INTERNAL  NODE  MOMENT  FLAGS 

3280  C 

3290  DO  180  JL=2 ,KR0W 

3300  NSTART=KC0L-NN0DE( JL)+1 

3310  DO  170  JN  =  NSTART  ,KC0L 

3320  MV(JL,JN)=1. 

3330  MH(JL,JN)=1. 

3340  170      CONTINUE 

3350  180   CONTINUE 

3360  C234567 

71 


3370    C 

33  80    C   CREATE  A  LINE  ARRAY  OF  NON-ZERO  VERTICAL  MOMENTS  AND  AN 

3390    C   ASSOCIATED  ARRAY  OF  MAXIMUM  PLASTIC  MOMENTS;  COUNT 

3400    C   TOP  TO  BOTTOM,  LEFT  TO  RIGHT 

3410    C 

3420  NLV=1 

3430  DO  2 10  J=2 ,KCOL 

3440  DO  200  1=1 ,KROW 

3450  IF(MV(I, J)  .LT. 1  )  GO  TO  200 

3460  MLIN(NLV)=MV(I  ,  J) 

3470  MMAX(NLV)=VERTMO( J) 

3480  NLV=NLV+1 

3490     200      CONTINUE 

3500     210   CONTINUE 

3510    C 

3520    C   CONTINUE  THE  LINE  ARRAY  WITH  NON-ZERO  HORIZ.  MOMENTS  AND  THEIR 

3530    C   ASSOCIATED  LIMITS. 

3540    C 

3550    C234567 

3560  NLH=NLV 

3570  DO  230  1=2 ,KROW 

3580  DO  220  J=l  ,KCOL 

3590  IF(MH( I, J)  .LT.  1  )  GO  TO  220 

3600  MLIN(NLH)=MH( I  ,  J) 

3610  MMAX(NLH)=HORZMO(I) 

3620  NLH=NLH+1 

3630     220      CONTINUE 

3640     230      CONTINUE 

3650    C 

3660    C 

3670    C   CALCULATE  TOTAL  NUMBER  OF  EQUATIONS  (NODAL  EQUILIBRIUM  PLUS 

3680    C   INEQUALITY  EQUATIONS) 

3690  NEQT0T=N0DEEQ+NLH-1 

3700    C 

3710    C 

3720    C   ESTABLISH  THE  C  MATRIX  (OBJECTIVE  FUNCTION  VECTOR) 

3730    C   THIS  IS  THE  FINAL  C  MATRIX 

3740    C 

3750  DO  240  1=1 ,NLH 

3760  C(I)=0.0 

3770     240   CONTINUE 

3780  C(NLH)=1.0 

3790    C 

3800    C 

3810    C   INITIALLIZE  THE  AA  MATRIX  (VARIABLE  COEFFICIENTS) 

3820    C   NOTE:  NO.  OF  ROWS  REQUIRED  BY  ZX3LP  IS  NEQTOT+2 

3830    C 

3840    C 
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3850 

3860  DO    250    I=l,NEQT0T+2 

3870  DO    245    J=1,NLH 

3880  AA(I,J)=0.0 

3890  245  CONTINUE 

3900  250       CONTINUE 

3910  C234567 

3920  C 

3930  C   INITIALIZE  THE  B  MATRIX  (RIGHT  HAND  SIDE  VECTOR) 

3940  C   FILL  IN  ANY  KNOWN  (INEQUALITY)  VALUES  FOR  B 

3950  C   IF  BOUNDARIES  ARE  FIXED,  THIS  IS  THE  FINAL  B  MATRIX 

3960  C 

3970  DO    255     I=l,NEQT0T+2 

3980  IF(I.GT.NLH-1 )     GO    TO    254 

3  990  B(I)=2  .*MMAX(I) 

4000  GO    TO    255 

4010  254  B(I)=0.0 

4020  255       CONTINUE 

4030  C 

4040  C   FILL  IN  AA  MATRIX  INEQUALITY  CONSTRAINT  VARIABLE  COEFFICIENTS 

4050  C 

4060  DO  265  1=1 ,NLH-1 

4070  AA(I,I)=1. 

4080     265   CONTINUE 

4090  C 

4100  C   LOAD  NODAL  FORCE  COEFFICIENTS  INTO  THE  AA  MATRIX 

4110  C 

4120  INDEX=NLH 

4130  DO       270     I=1,N0DEEQ 

4140  AA(INDEX,NLH)=    -1.0 

4150  INDEX=INDEX+1 

4160     270   CONTINUE 

4170  C 

4180  C   CREATE  A  LINE  ARRAY  OF  HORIZ.  LENGTHS,  AND  ANOTHER  ARRAY 

4190  C   OF  VERTICAL  LENGTHS 

4200  C 

42  10  C  HORIZONTAL    ARRAY: 

4220  C 

4230  NNN=1 

4240  DO    2  80    J=2  ,KR0W 

4250  NSTART=KCOL-NNODE( J)+l 

4260  DO    275    L=NSTART , KCOL 

42  7  0  HORZL(NNN)  =  XCOORD( J  ,L)-XCOORD( J  ,L-1 ) 

42  80  NNN=NNN+1 
4290            275                CONTINUE 
4300             280       CONTINUE 

43  10  C 

4320  C   VERTICAL  ARRAY: 
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4330  NNN=1 

4340  DO    290    J=2 ,KCOL 

4350  NSTART=KR0W-NVERT(J)+1 

4360  DO    285     I=NSTART ,KR0W 

43  7  0  VERTL(NNN)=YC00RD(I-1  , J ) -YCOORD ( I , J  ) 

43  80  NNN=NNN+1 

4390  285               CONTINUE 

4400  290       CONTINUE 

4410  C 

4420  C   CREATE  AN  ARRAY  ( KHORZ ( N ) / /KVERT ( N ) )  WHICH  ASSOCIATES  THE 

4430  C   SEQUENTIAL  ( HORI ZONTAL / / VERT ICAL)  NODE  NUMBER  TO  EACH  INTERNAL 

4440  C   NODE  NUMBER.  INTERNAL  NODE  NUMBERING  LEFT-RIGHT  ,  TOP-BOTTOM 

4450  C 

4460  NINT=1 

4470  DO  300  1=2 ,KR0W 

4480  NSTART=KC0L-NN0DE(I)+1 

4490  DO  295  J=NSTART , KCOL 

^500  IF(IBC.EQ.l)  KHORZ(NINT)=NINT 

^510  IF(IBC.EQ.2)  KHORZ ( NINT )  =  NINT+( I - 1  ) 

4  52  0  C2 34567890123456789012345678901234567890123456789012345678901234567890 

4530  IF(J.EQ.2)  K0UNT=0 

4540  IF(J.GT.2)  GO  TO  296 

^550  298           KVERT(NINT)=K0UNT+I+IBC-1-(KR0W-NVERT(J) ) 

4560  GO  TO  299 

4570  296           KOUNT=0 

4580  DO  297  K=2 , J-1 

^590  K0UNT=K0UNT+NVERT(K)+IBC-1 

4600  297           CONTINUE 

4610  GO  TO  298 

4620  299           NINT=NINT+1 

4630  295      CONTINUE 

4640  300   CONTINUE 

4650  C 

4660  C   CREATE  A  LINEAR  ARRAY  OF  HORIZ.  AND  VERT.  PLASTIC  MOMENTS 

4670  C   ASSOCIATED  WITH  EACH  INTERNAL  NODE 

4680  C 

4690  NINT=1 

4700  DO  3 10  1=2 ,KR0W 

4710  NSTART=KC0L-NN0DE(I )+l 

^720  DO  305  J=NSTART,KC0L 

4730  MPV(NINT)=VERTMO(J) 

4740  MPH(NINT)=H0RZM0(I) 

4750  NINT=NINT+1 

4760  305      CONTINUE 

4770  310   CONTINUE 

4780  C 

4790  C   CREATE  AN  ARRAY  ASSOCIATING  VERTICAL  LENGTH  INDICES  WITH 

4800  C   INTERNAL  NODES 
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4810  C 

4820  NINT=1 

^830  DO  330  1=2 , KROW 

^840  NSTART=KCOL-NNODE(I )+l 

^850  DO  325  J=NSTART , KCOL 

^860  IF(J.EQ.2)     KOUNT=0 

^870  IF( J.GT.2  )    GO    TO    326 

^880  328                           NVLEN(NINT)=KOUNT+I-(KROW-NVERT(J)) 

4890  GO    TO    329 

4900  326                          KOUNT=0 

^910  DO  327  K=2 , J-1 

^920  K0UNT=K0UNT+NVERT(K) 

4930  327           CONTINUE 

4940  GO  TO  328 

4950  329          NINT=NINT+1 

4960  325        CONTINUE 

4970  330   CONTINUE 

4980  C 

4990  C    CREATE  AN  ARRAY  OF  FLAGS  ASSOCIATED  WITH  EACH  INTERNAL  NODE 

5000  C    TO  INDICATE  THE  FOLLOWING  FLAG  CONDITION: 

5010  C        (1)  NODE  IS  NEAR  VERTICAL  BOUNDARY 

5020  C        (2)  NODE  IS  NEAR  HORIZONTAL  BOUNDARY 

5030  C        (3)  NODE  IS  NEAR  BOTH  BOUNDARIES 

5040  C        (4)  NODE  IS  A  SYMMETRY   NODE  ABOUT  THE  HORIZ.  DIAMETER 

5050  C        (5)  NODE  IS  A   SYMMETRY  NODE  ABOUT  THE  VERT.  DIAMETER 

5060  C        (6)  NODE  IS  NEAR  BOUNDARY  AND  NEAR  HORIZONTAL  AXIS  OF  SYMMETRY 

5070  C        (7)  NODE  IS  NEAR  BOUNDARY  AND  NEAR  VERTICAL  AXIS  OF  SYMMETRY 

5080  C 

5090  NINT=1 

5100  DO  320  1=2 ,KROW 

5110  NSTART=KCOL-NNODE(I)+I 

5120  DO  315  J=NSTART,KCOL                                             , 

5130  IF(HORZL(NINT)  .LT.XDIST)  I  FLAG ( NINT ) =2 

5140  IF(VERTL(NVLEN(NINT)) .LT.YDIST)  IFLAG ( NINT ) = 1 

^^5°  IF ( VERTL ( NVLEN ( NINT ) ).LT. YD  I  ST. AND. HORZL( NINT). LT.XDIST) 

5160  C23 4567 890 123 4567 890 123 4567 890 123 4567890 123 4567 890 123 4567 890 123 4567 890 

5170  1   IFLAG(NINT)=3 

5180  IF(J.EQ.KCOL)  IFLAG ( NINT )= 5 

5190  IF(J.EQ.KC0L.AND.I.EQ.2 )  IFLAG ( NINT )= 7 

5200  IF(I.EQ.KROW)  IFLAG ( NINT ) =4 

5210  IF(I.EQ.KROW.AND.J.EQ.NSTART)  I  FLAG ( NINT )  =  6 

5220  NINT=NINT+1 

5230  315      CONTINUE 

5240  320   CONTINUE 

5250  DO  1000  I=1,N0DEEQ 

5260  WRITE(6,1001)  HORZ L ( I ) , NVLEN ( I ) , VERTL ( NVLEN ( I ) ) 

5270  1000  CONTINUE 

5280  1001  F0RMAT(1X,F8.3 ,I8,F8.3 ) 
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5290  C234567 

5300  C 

53  10  C  LOAD  NODAL  EQUATION  COEFFICIENTS  INTO  AA  MATRIX 

5320  C 

5330  INDEX=NLH-1 

5340  DO    340    I=1,N0DEEQ 

5350  IF(IBC.EQ.  1)  GO  TO  341 

5360  KHVAL=NLV-1+KH0RZ(I) 

5370  KVVAL=KVERT(I) 

5380  IF(I.EQ.NODEEQ)  GO  TO  342 

5390  AA(INDEX+I,KHVAL-1)=  -l./HORZL(I) 

5^00  IF(IFLAG(I) .EQ.5.AND.I0DDV.EQ.1)  GO  TO  343 

5^10  IF(IFLAG(I) .EQ.5.AND.IODDV.EQ.0)  GO  TO  344 

l^l^^  IF(IBC.EQ.2.AND.IFLAG(I).EQ.7.AND.IODDV.EQ.l)  GO  TO  3  43 

l,]l  IF(IBC.EQ.2.AND.IFLAG(I).EQ.7.AND.IODDV.EQ.O)  GO  TO  344 

^^t°  AA(INDEX+I,KHVAL)=  l./HORZL(I)+  1 . / HORZL ( I+l ) 

5f50  AA(INDEX+I,KHVAL+1)=  - 1 ./ HORZL ( I+l ) 

5460  GO  TO  345 

5^70  343       AA(INDEX+I,KHVAL-1)=  -2 ./HORZL(I) 

5^80  AA(INDEX+I,KHVAL)=  2./H0RZL(I) 

5490  GO  TO  345 

5500  344      AA(INDEX+I,KHVAL-1)=  -l./HORZL(I) 

5510  AA(INDEX+I,KHVAL)=  l./HORZL(I) 

5520  KVVAL=KVERT(I) 

5530  345       IF(I  .EQ.NODEEQ)  GO  TO  342 

5540  AA(INDEX+I,KVVAL-1)=  - 1 . / VERTL ( NVLEN ( I ) ) 

^^50  IF(IFLAG(I) .EQ.4.AND.I0DDH.EQ.  1)  GO  TO  346 

^5^°  IF(IFLAG(I) .EQ.4.AND.I0DDH.EQ.0)  GO  TO  347 

VCJ^  IF(IBC.EQ.2.AND.IFLAG(I).EQ.6.AND.IODDH.EQ.l)  GO  TO  346 

iiTi  IF(IBC.EQ.2.AND.IFLAG(I).EQ.6.AND.IODDH.EQ.O)  GO  TO  347 

llll  AA(INDEX+I,KVVAL)=  1 . /VERTL ( NVLEN ( I ) )+  1 ./ VERTL ( NVLEN ( I ) +1) 

5600  C234567890123 4567890  123 4567890123 4567890123 4567890123 4567890 123 4567890 

5610  AA(INDEX+I,KVVAL+1)=  - 1 ./ VERTL ( NVLEN ( I )+l ) 

5620  GO  TO  340 

^"0  346       AA(INDEX+I,KVVAL-1)=  -2  .  /  VERTL  (  N  VLEN  (  I  )  ) 

^?^°  AA(INDEX+I,KVVAL)=  2 . / VERTL ( NVLEN ( I ) ) 

5650  GO  TO  340 

5660  347      AA(INDEX+I,KVVAL-1)=  - 1 ./ VERTL ( NVLEN ( I ) ) 

^^'^  AA(INDEX+I,KVVAL)=  1 ./ VERTL ( NVLEN ( I ) ) 

5680  GO  TO  340 

5690  342       IFdODDV.EQ.l  )G0  TO  348 

5^00  AA(INDEX+I,KHVAL)=  l./HORZL(I) 

11^'^  AA(INDEX+I,KHVAL-1)=  -l./HORZL(I) 

5720  GO  TO  349 

5*^30  348      AA(INDEX+I,KHVAL)=  2./H0RZL(I) 

^-^^^  AA(INDEX+I,KHVAL-1)=  -2./H0RZL(I) 

5750  349       IF(IODDH.EQ.O)  GO  TO  350 

5760  AA(INDEX+I,KVVAL)=  1 ./ VERTL ( NVLEN ( I ) ) 
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5770  AA(INDEX+I ,KVVAL-1 )=  - 1 . / VERTL ( NVLEN ( I ) ) 

5780  GO  TO  340 

5790  350      AA(INDEX+I ,KVVAL)=  2 . / VERTL ( N VLEN ( I ) ) 

5800  AA(INDEX+I,KVVAL-1 )=  -2 . / VERTL ( NVLEN ( I ) ) 

5810  340      CONTINUE 

5820  IF(KKK.EQ.l)  GO  TO  341 

5830  WRITE(18,3000) 

5840  WRITE( 18,3001 )  ( ( AA( I , J )  , J=l  , NLH )  ,  1  =  1  ,  NEQTOT+2 ) 

5850  3000     F0RMAT( IX, 'THE  AA  MATRIX  IS:') 

5860  3001     F0RMAT(1X, 13 (2X,F7.4) ) 

5870  341   PRINT  *,  '   SIMPLE  SUPPORT  NOT  FULLY  IMPLEMENTED  YET' 

5880  IF(IBC.EQ.l)  STOP 

5890  M1=NLH-1 

5900  M2=N0DEEQ 

5910  N=NLH 

5920  IA=M2+MH-2 

5930  IRWl=IA*IA+3*Ml+2*M2+4 

5940  IWl=2*M2+3*Ml+4 

5950  MN=MAX0(N,M1+M2 ) 

5  9  60  CALL    DIMEN(A,AA,B , C , N ,M 1  ,M2  , I A , IRW 1 , IW 1  ,  PSOL , DSOL , IW , RW ,MN) 

5970  C 

5980  C  CONVERT  ZX3LP  SOLUTION  VECTOR  PSOL(I)  TO  THE  CORRESPONDING 

5990  C  NORMALIZED  MOMENTS  AT  THE  NODES,  AND  ALSO  DETERMINE  THE 

6000  C  APPLIED  NODAL  FORCE  AND  EQUIVALENT  UNIFORM  PRESSURE 

6010  C 

6020  DO  400  I=1,NLH-1 

603  0  MLIN(I)=(PSOL(I)-MMAX( I) )/MMAX(I) 

6040  400   CONTINUE 

6050  IF(HFLAG.EQ. 1 )  HFAC=0 

6060  IF(HFLAG.EQ.O)  HFAC=1 

6070  IF(VFLAG.EQ. 1 )  IFAC=0 

6080  IF(VFLAG.EQ.O)  IFAC=1 

6090  FN0DE=PS0L(NLH) 

6100  FT0T=FN0DE*(4 .*N0DEEQ-HFAC*2 . *NN0DE ( KROW ) -IF AC*2 . *NVERT ( KCOL ) 

6110  1  +HFAC*IFAC) 

6120  PRESS=FT0T/(3 . 14159265*R*R) 

6130  WRITE(18,450)  PRESS 

6140  WRITE( 18 ,460) 

6150  WRITE( 18,470) 

6160  DO  500  I=1,NLV-1 

6170  WRITE( 1 8, 480)  I , MLIN ( I ) , MMAX( I ) 

6180  500   CONTINUE 

6190  WRITE( 18,490) 

6200  WRITE(18,485) 

6210  DO  600  I=NLV,NLH-1 

6220  WRITE( 18,480)  I-NLV+1 ,MLIN ( I ) ,MMAX( I ) 

6230  600   CONTINUE 

6240  450   FORMAT( IX, 'THE  COLLAPSE  PRESSURE  IS:',F10.2,'  PSI') 
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6250  460   FORMATdX,  'THE  NORMALIZED  VERTICAL  MOMENTS  ARE: 

6260  1(T0P-B0TT0M, LEFT-RIGHT) ' ) 

6270  470   FORMATdX,  'VERT  MOMENT  NO  .',  5X  ,' NORMAL  IZ  ED  VALUE', 5X, 

6280  I'MAX  PLASTIC  MOMENT') 

6290  480   F0RMAT(8X,I3 ,9X,F17.4,5X,E15.4) 

6300  485   FORMATdX, 'HORZ  MOMENT  NO  .',  5X  ,' NORMAL  IZED  VALUE', 5X 

6310  I'MAX  PLASTIC  MOMENT') 

6320  490   FORMATdX, 'THE  NORMALIZED  HORIZONTAL  MOMENTS  ARE: 

6330  1 (LEFT-RIGHT, TOP-BOTTOM) ' ) 

63  40  C23 4567 890 123 4567 890 123 4567 890 123 4567 890 123 456 7 890 123 456 7 890 123 4567  891 

6350  STOP 

6360  END 

63  70 

6380 

6390  C234567 

6400 

6410  SUBROUTINE  DIMEN ( A , AA , B , C , N , Ml ,M2 , I A , IRW 1 , IW 1 , PSOL 

6A20  1  DSOL, IW,RW,MN) 

6^3  0  DIMENSION  AAd2  0,80),A(IA,N),B(IA),C(N),PSOL(MN) 

6440  DIMENSION  DSOL ( lA ), IW( IW 1 ), RW ( IRW 1 ) 

6450  DO  12  J=1,N 

6460  DO  11  I=l,IA-2 

6470  A(I, J)=AA(I ,J) 

6480  11     CONTINUE 

6490  12  CONTINUE 

11^^  CALL  ZX3LP(A,IA,B,C,N,M1,M2 , S , PSOL , DSOL , RW , IW , lER) 

6510  RETURN 

6520  END 

6530  C 

65^0  C       THIS  SUBPROGRAM  CALCULATES  PLASTIC  MOMENT  OF  EQUIVALENT 

6550  C      STIFFENER.  WHEN  CALLING  THIS  FUNCTION ," SPACE"  IS  THE 

6560  C       STIFFENER  SPACING  BETWEEN  ADJACENT  STIFFENERS. 
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100 

C 

no 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

22  0 

230 

240 

250 

260 

C 

270 

c 

2  80 

c 

290 

c 

300 

c 

310 

320 

c 

330 

340 

350 

360 

FUNCTION  MP(DW,SPACE ,T2 ,B1 ,T1 ,SIGMAY,T3 ) 

REAL  MP 

R1=(DW*T3 )/(SPACE*T2 ) 

R2=(B1*T1 ) /(SPACE*T2 ) 

R3=T2/DW 

R4=T1/DW 

IF(R1+R2  .LT  .0  .99999)  GO  TO  5 

Z  =  SPACE*DW*T2*(Rl*(RH-2  .  *R3 +2  .  ) +2  .  *R1  *R2  *  (  R4  +  1  .  )  -  (R2  -  1  .  )  **2  ) 
1/(4. *R1) 

GO  TO  6 

Z=SPACE*DW*T2*(2 . *Rl+4 . *R2 +R3 * ( 1 .+2 . *R 1+2 . *R2 )-R3*(Rl+R2 )**2 
1+2 .*R2*R4)/4. 

MP=Z*SIGMAY 

RETURN 

END 


THIS  SUBPROGRAM  CREATES  BLANK  LINES  ON  SCREEN 


C234567 

FUNCTION    SKIP(N) 
DO     10    1=1, N 

SKIP=I 

PRINT    *,     '        * 
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no 

120 


10 


CONTINUE 

RETURN 

END 
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INTENTIONALLY  BLANK 
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INTENTIONALLY  BLANK 
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APPENDIX  B    ZX3LP 


ZX3LP  LINEAR  PROGRAM. 


This  section,  taken  from  reference  {15},  is  intended  to 
provide  a  brief  explanation  of  the  capabilities  and  uses  of 
the  linear  programming  subroutine  ZX3LP. 

IMSL  ROUTINE  NAME   -  ZX3LP 

PURPOSE  -  SOLVE  THE  LINEAR  PROGRAMMING  PROBLEM 

VIA  THE  REVISED  SIMPLEX  ALGORITHM, - 
EASY  TO  USE  VERSION 

USAGE  -  CALL  ZX3LP  ( A , I A , B ,C ,N,M1 ,M2 , S , PSOL, 

DSOL,RW,IW,IER) 

ARGUMENTS      A     -  MATRIX  OF  DIMENSION  M1+M2+2  BY  N 

CONTAINING  THE  COEFFICIENTS  OF  THE  Ml 
INEQUALITY  CONSTRAINTS  IN  THE  FIRST  Ml 
ROWS  FOLLOWED  BY  THE  COEFFICIENTS  OF 
THE  M2  EQUALITY  CONSTRAINTS.  (INPUT) 
THE  LAST  TWO  ROWS  OF  A  ARE  USED  ONLY 
AS  WORKING  STORAGE. 
I A    -  ROW  DIMENSION  OF  MATRIX  A  EXACTLY  AS 

SPECIFIED  IN  THE  DIMENSION  STATEMENT 
IN  THE  CALLING  PROGRAM.  (INPUT) 
TWO  ROWS  OF  A  ARE  REQUIRED  FOR  WORKING 
STORAGE,  AND  THEREFORE,  I A  MUST 
NOT  BE  LESS  THAN  M1+M2+2. 
B     -  VECTOR  OF  LENGTH  M1+M2+2  CONTAINING  THE 
RIGHT  HAND  SIDES  OF  THE  INEQUALITY 
CONSTRAINTS  IN  ITS  FIRST  Ml  LOCATIONS 
FOLLOWED  BY  THE  M2  RIGHT  HAND  SIDES  OF 
THE  EQUALITY  CONSTRAINTS.  (INPUT) 
THE  LAST  TWO  ELEMENTS  OF  B  ARE  USED 
AS  WORKING  STORAGE. 
C     -  VECTOR  OF  LENGTH  N  CONTAINING  THE 

COEFFICIENTS  OF  THE  OBJECTIVE  FUNCTION 
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(INPUT) 
N     -  NUMBER  OF  UNKNOWNS  IN  THE  MODEL.  (INPUT) 
Ml    -  NUMBER  OF  INEQUALITY  CONSTRAINTS .( INPUT) 
M2    -  NUMBER  OF  EQUALITY  CONSTRAINTS.   (INPUT) 
S     -  VALUE  OF  THE  OBJECTIVE  FUNCTION . (OUTPUT) 
PSOL  -  VECTOR  OF  LENGTH  N  CONTAINING  THE  PRIMAL 
SOLUTION.   (OUTPUT)  PSOL  IS  ALSO  USED 
AS  WORK  STORAGE  AND  THEREFORE  MUST 
HAVE  LENGTH  AT  LEAST  MAX(N ,M1+M2 ) . 
DSOL  -  VECTOR  OF  LENGTH  M1+M2+2  CONTAINING  THE 
DUAL  SOLUTION.   (OUTPUT) 

THAT  IS,  DSOL(l) ,DS0L(M1+M2) 

CONTAIN  THE  SOLUTION  TO  THE 
PROBLEM  MIN  BT*Y  SUBJECT  TO  AT*Y  IS 
GREATER  THAN  OR  EQUAL  TO  C  AND  Y 
GREATER  THAN  OR  EQUAL  TO  0  WHERE 
AT=A-TRANSPOSE  AND  BT=B-TRANSPOSE . 
WHEN  THE  PRIMAL  PROBLEM  HAS 
EQUALITY  CONSTRAINTS,  THE 
CORRESPONDING  COMPONENTS  OF  THE 
DUAL  SOLUTION  ARE  UNCONSTRAINED. 
DS0L(M1+M2+1)  AND  DSOL(Ml+M2+2 ) 
ARE  USED  AS  WORKING  STORAGE. 
RW    -  WORK  VECTOR  OF  LENGTH 

(Ml+M2+2)*(Ml+M2+2)+  3*Ml+2*M2+4. 
IW    -  WORK  VECTOR  OF  LENGTH  2*M2+3*Ml+4. 
lER   -  ERROR  INDICATOR.  (OUTPUT) 
TERMINAL  ERROR 
lER  =  130  INDICATES  THAT  I A  IS  LESS 

THAN  M1+M2+2. 
lER  =  131  INDICATES  THAT  THE  COST 
CRITERION  HAS  UNBOUNDED  VALUES. 
lER  =  132  INDICATES  THAT  THE  MAXIMUM 
NUMBER  OF  ITERATIONS  WAS  REACHED  IN 
ZXOLP. 
lER  =  133  INDICATES  THAT  NO  FEASIBLE 
SOLUTION  EXISTS. 
WARNING  (WITH  FIX) 
lER  =  7  0  INDICATES  THAT  SOME 

ARTIFICIAL  VARIABLES  REMAINED  IN  THE 
SOLUTION  BASIS  AT  A  ZERO  LEVEL  AFTER 
PHASE  1.  THIS  CONDITION  CAN  BE  CAUSED 
BY  HAVING  REDUNDANT  CONSTRAINTS. 
NEVERTHELESS,  A  SOLUTION  IS  COMPUTED 
AND  RETURNED  IN  PSOL  AND  DSOL. 
PRECIS I ON /HARDWARE  -  SINGLE  AND  DOUBLE/H32 

-  SINGLE/H36,H48,H60 
REQD.  IMSL  ROUTINES-  UERTST,UGETIO, ZXOLP 
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NOTATION  -  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE 
UHELP. 
REMARKS   ZX3LP  IS  INTENDED  TO  VERY  EASY  TO  USE.   THEREFORE. 

THE  CALLING  PARAMETERS  ARE  VERY  SIMPLE  AND  THERE  ARE 
NO  OPTIONS.   ZXOLP  CAN  BE  USED  DIRECTLY  IN  SITUATIONS 
THAT  REQUIRE  MORE  FLEXIBILITY  THAN  IS  PROVIDED  BY 
ZX3LP. 


ALGORITHM 

To  solve  the  linear  programming  problem, 
subject  to 


maximize   c-^PSOL^    +    ...+   Cj^PSOL^^    =    SP 


il  PSOLi    +    ...    +   aiNPSOL<b^  i  =  l,...,Mi 

a^iPSOLi    +    ...    +    aiNPSOLN  =  bi  i=M   +1,  ...,M 

PSOLj>0  j=l,...,N 

where   M   =   Mi    +   Mo. 

The  dual  linear  programming  problem  is, 

minimize  b-[^DSOL;L  +•.."*■  bj^DSOLp^  =  SD 
subject  to 

a,-DSOL,  +  ...  +  aMjDSOLM^c^       j=l,...,N 
DSOLi^O        i=l,...,Mi 
DSOL^  unrestricted  in  sign  when    i=Mi+l,...,M 

ZX3LP  computes  the  solution  to  the  primal  problem,  PSOL, 
the  solution  to  the  dual  problem,  DSOL,  and  the  values  of 
the  objective  function  S=SP=SD. 

ZX3LP  calls  ZXOLP  which  solves  the  linear  programming 
problem  by  the  revised  simplex  method. 

See  reference  {14}  for  additional  detail. 


PROGRAMMING  NOTES 


ZX3LP   is   intended  to  very  easy  to  use.   Therefore,  the 
calling  parameters  are  very  simple  and  there  are  no  options. 
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ZXOLP  can   be  used  directly  in  situations  that  require  more 
flexibility  than  is  provided  by  ZX3LP. 


EXAMPLE 


Suppose  we  want  to  maximize   x 

subject  to  X-, 


-XI  -  X2 

X;L>0  ,X2 


+  3x  =  S 

<  1 

<  1 

<  1.5 

<  -.5 

>  0 


X2 


We  can  proceed  as  follows: 

INTEGER    IA,N,M1,M2,IW(16) ,IER 

REAL       A(6,2),B(6),C(2),RW(52) , PSOL ( 4 ) ,DS0L(6) ,S 

N   =   2 

Ml  =   4 

M2  =   0 

lA  =   6 


1.  0. 

A   =     0.  1. 

1.  1. 

•1.  -1. 

'l.  ' 
1. 
1.5 
-.5 


C   =     (1.   3.) 

CALL  ZX3LP  (A,IA,B,C,N,M1,M2,S,PS0L,DS0L,RW,IW,IER) 


END 

Output 

lER  =  0 

S    =  3.5 

PSOL  =  (  .  5 

1.) 

DSOL  =  (0. 

2.   1. 

0.) 
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INPUT  INFORMATION 

NO.  OF  VERTICAL  STIFFENERS:      10 

NO.  OF  HORIZONTAL  STIFFENERS:      10 

PLATE  DIAMETER:      lO.OO(FEET) 

VERTICAL  BEAMS 
MATERIAL  YIELD  STRESS:         50000.00  PSI 
PLATE  THICKNESS:        2.000  INCHES 
COMMON  WEB  THICKNESS:        0.000  INCHES 
COMMON  WEB  DEPTH      :        1.000  INCHES 
COMMON  BOTTOM  FLANGE  THICKNESS:        0.000  INCHES 
COMMON  BOTTOM  FLANGE  WIDTH      :        0.000  INCHES 


HORIZONTAL  BEAMS 

COMMON  WEB  THICKNESS:  0.000  INCHES 

COMMON  WEB    DEPTH   :  1.000  INCHES 
COMMON  BOTT .FLG. THICKNESS :        0.000  INCHES 
COMMON  BOTT. FLG.   WIDTH   :        0.000  INCHES 


OUTPUT  INFORMATION 

UNSUPP.  PLATE  SECTION  COLLAPSE  LOAD :  1  4 72 2  .  22 2 7 
THE  COLLAPSE  PRESSURE  IS:      154.70  PSI 

THE  NORMALIZED  VERTICAL  MOMENTS  ARE :( TOP-BOTTOM , LEFT-RIGHT ) 

VERT  MOMENT  NO.  NORMALIZED  VALUE       MAX  PLASTIC  MOMENT 

1  -1.0000  0.5455E+06 

2  0.0752  0.5455E+06 

3  1.0000  0.5455E+06 

4  1.0000  0.5455E+06 

5  -1.0000  0.5455E+06 

6  -0.0847  0.5455E+06 

7  1.0000  0.5455E+06 

8  1.0000  0.5455E+06 

9  1.0000  0.5455E+06 

10  -1.0000  0.5455E+06 

11  -1.0000  0.5455E+06 

12  1.0000  0.5455E+06 

13  1.0000  0.5455E+06 

14  0.8795  0.5455E+06 

15  1.0000  0.5455E+06 

16  -1.0000  0.5455E+06 

17  -0.9338  0.5455E+06 

18  -1.0000  0.5455E+06 

19  -0.1566  0.5455E+06 

20  1.0000  0.5455E+06 
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21 

-0.6024 

0.5455E+06 

22 

-1 .0000 

0.5455E+06 

23 

-0.3511 

0.5455E+06 

24 

-0.0847 

0.5455E+06 

25 

0.0571 

0.5455E+06 

26 

-0.1988 

0.5455E+06 

27 

1 .0000 

0.5455E+06 

THE  NORMALIZED 

HORIZONTAL  MOMENTS  ARE: 

: (LEFT-RIGHT, TOP-BOTTOM) 

HORZ  MOMENT  NO. 

NORMALIZED  VALUE 

MAX  PLASTIC  MOMENT 

1 

-1 .0000 

0.5455E+06 

2 

0.7521 

0.5455E+06 

3 

1.0000 

0.5455E+06 

4 

1  .0000 

0.5455E+06 

5 

-1  .0000 

0.5455E+06 

6 

-0.  83  10 

0.5455E+06 

7 

-0.8532 

0.5455E+06 

8 

0.7270 

0.5455E+06 

9 

1 .0000 

0.5455E+06 

10 

-1  .0000 

0.5455E+06 

11 

-1.0000 

0.5455E+06 

12 

-0.6990 

0.5455E+06 

13 

0.2891 

0. 5455E+06 

14 

1  .0000 

0.5455E+06 

15 

1  .0000 

0. 5455E+06 

16 

-1  .0000 

0.5455E+06 

17 

-1 .0000 

0.5455E+06 

18 

-0,4728 

0.5455E+06 

19 

-0.3433 

0.5455E+06 

20 

-0.8523 

0.5455E+06 

21 

1  .0000 

0.5455E+06 

22 

-1  .0000 

0.5455E+06 

23 

0.0480 

0.5455E+06 

24 

0.7228 

0.5455E+06 

25 

1 .0000 

0.5455E+06 

26 

1  .0000 

0.5455E+06 

27 

-I  .0000 

0.5455E+06 
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INPUT  INFORMATION 

NO.  OF  VERTICAL  STIFFENERS:       4 

NO.  OF  HORIZONTAL  STIFFENERS:       5 

PLATE  DIAMETER:      12.00(FEET) 

VERTICAL  BEAMS 
MATERIAL  YIELD  STRESS:         50000.00  PSI 
PLATE  THICKNESS:        1.000  INCHES 
COMMON  WEB  THICKNESS:        0.800  INCHES 
COMMON  WEB  DEPTH     :        5.000  INCHES 
COMMON  BOTTOM  FLANGE  THICKNESS:        1.000  INCHES 
COMMON  BOTTOM  FLANGE  WIDTH     :        4.000  INCHES 


HORIZONTAL  BEAMS 

COMMON  WEB  THICKNESS:  0.750  INCHES 

COMMON  WEB    DEPTH   :  6.000  INCHES 
COMMON  BOTT. FLG. THICKNESS :        1.000  INCHES 
COMMON  BOTT. FLG.   WIDTH   :        4.000  INCHES 


OUTPUT  INFORMATION 

UNSUPP.  PLATE  SECTION  COLLAPSE  LOAD:   588.3488 
THE  COLLAPSE  PRESSURE  IS:      142.43  PSI 

THE  NORMALIZED  VERTICAL  MOMENTS  ARE :( TOP-BOTTOM , LEFT-RIGHT ) 

MAX  PLASTIC  MOMENT 
0.2132E+07 
0.2132E+07 
0.2132E+07 
0.2 132E+07 
0.2132E+07 
0.2  132E+07 
0.2132E+07 
0.2132E+07 
(LEFT- RIGHT , TOP -BOTTOM) 
MAX  PLASTIC  MOMENT 
0.2450E+07 
0.2450E+07 
0.2450E+07 
0.2450E+07 
0.2450E+07 
0.2450E+07 
0.2450E+07 
0.2450E+07 
0  .2450E+07 


VERT  MOMENT 

NO. 

NORMALIZED  VALUE 

1 

-1.0000 

2 

-1 .0000 

3 

-0.4645 

4 

1 .0000 

5 

-1  .0000 

6 

0. 1584 

7 

1  .0000 

8 

0.53  62 

THE  NORMALIZED 

HORIZONTAL 

MOMENTS  ARE 

HORZ  MOMENT 

NO. 

NORMALIZED  VALUE 

1 

-1  .0000 

2 

0.0455 

3 

1  .0000 

4 

-1 .0000 

5 

1  .0000 

6 

1 .0000 

7 

-1 .0000 

8 

-0.8478 

9 

1  .0000 

90 
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